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ABSTRACT 

An  application  of  two-dimensional  Fourier  Transform  correlation 
for  the  attitude  measurement  of  three-dimensional  objects  is  considered, 
Included  are  geometric  optical  equations  and  geometric  translation  and 
rotation  equations  necessary  for  the  image  error  calculations.  Mini- 
mum attitude  resolution  equations  are  developed  to  estimate  image 
quantization  necessary  for  a  given  measurement  error.  These  techniques 
are  then  applied  as  an  additional  measurement  to  the  radar  tracking 
problem. 
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I.   INTRODUCTION 

The  development  of  sophisticated  naval  weapons  has  resulted  in  the 
requirement  for  high  accuracy  measurement  tests  and  discrete  analysis  of 
complex  system  parameters.  Typical  of  this  precise  test  and  evaluation 
requirement  is  attitude  measurement  of  three-dimensional  airborne  objects 
by  processing  multiple  two-dimensional  images. 

The  most  common  method  used  to  estimate  attitude  is  by  making 
absolute  or  precisely  relative  three-dimensional  geometric  trajectory 
measurements  of  target  object  extreme  structural  features  (nose,  tail, 
right  wing  tip,  left  wing  tip,  etc.)  from  camera  images,  By  appropriate 
camera,  range  coordinate  system,  and  target  coordinate  system  trans- 
formations these  measurements  will  yield  an  attitude  measurement. 

A  more  direct  method  of  attitude  measurement  involves  the  correlation, 
or  comparison,  of  an  image  of  an  object  at  an  unknown  attitude  with  a 
second  image  of  a  similar  object  at  a  known  attitude.  At  the  Naval 
Weapons  Center,  an  operator  controlled  film  and  television  correlation 
assessment  technique  (FATCAT)  facility  use  this  more  direct  technique 
for  attitude  estimation  [1]. 

This  thesis  presents  an  enhancement  of  the  correlation/comparison 
method  of  attitude  measurement  by  fully  automating  the  technique  incorp- 
orating a  pattern  recognition  method  of  the  directly  digitized  image. 
Image  theory,  resolution,  and  filtering  leading  to  the  two  dimensional 
matched  filter  attitude  measurement  concept  is  discussed  and  this  concept 
is  applied  as  an  additional  measurement  to  a  radar  tracking  problem. 
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II.   IMAGE  THEORY 

Three-dimensional  to  two-dimensional  projections  are  obtained  by 
viewing  the  real  world  as  recorded  by  a  camera  or  other  imaging  device. 
Image  theory  defines  projection  as  the  transformation  of  each  point  in 
the  field  of  view,  through  a  focal  point,  and  onto  a  plane.  This  is 
also  equivalent  to  locating  the  focal  plane  in  front  of  the  focal  point 
and  changing  the  sign  of  the  projected  points  [2].  These  geometrical 
relationships  are  illustrated  in  Fig.  1.  The  transformation  of  point 
X,  y,  z  in  the  X,  Y,  Z  coordinate  system  of  Fig.  1  is  given  by  the 
following  equations: 


u  = 


y  +  f 


(1) 


v=  z  (2) 

y  +  f 


Where  u  =  projected  x  dimension 
V  =  projected  z  dimension 

If  f  <<  y,  then  the  following  approximations  hold; 

fx 


u  = 


=  kx  (3) 


fz 
V  =  -y-  ky  (4) 


Therefore,  for  a  given  focal  length  f,  and  target  range  y,  a  constant 
image  scaling  factor  k,  is  calculated. 
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Figure  1 
Camera  Geometry 


A  second  coordinate  system  X',  Y',  Z'  can  be  used  to  describe  a 
three-dimensional  object  in  the  viewing  space.  If  the  origins  of  these 
two  coordinate  systems  are  coincident,  and  if  all  the  axes  are  parallel, 
then  the  transformation  of  points  from  one  coordinate  system  to  the 
other  would  be  simple. 


(5) 


x'  X 

y'  =  y 

z'  z 

-  i  L  . 
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If  the  axes  of  the  two  coordinate  systems  are  rotated  by  an  elevation 
angle  9  around  the  camera  Y  axis,  by  an  azimuth  angle  ^  around  the  camera 
Z  axis,  and  by  roll  angle  <t>   around  the  camera  X  axis,  the  transformation 
equation  would  be  [3],  [4]: 


x' 

COS  9 

0 

sin  8 

cos  ^ 

-sin  \i) 

0 

y' 

= 

0 

1 

0 

sinij^ 

COSip 

0 

z' 

-sine 

0 

cose 

0 

0 

1 

1    0      0 
0  COS  <p      -sin  (|) 
0  sin  4)   cos  4) 


(6) 


If  in  addition  there  was  a  translation  of  (AjB,C)  of  the  origin  of  the 
object  coordinate  system  relative  to  the  origin  of  the  camera  coordinate 
system,  the  transformation  equation  would  be: 


'x'' 

' 

' 

x" 

y' 

= 

e 

^ 

4) 

y 

z' 

z 

where 


cose    0 

sin  e 

8 

= 

0      1 

0 

- 

-sine    0 

cos  6 

- 

cos  ^ 

■sinijj      0 

^ 

= 

sin  tj; 

cos  ^      0 

- 

0 

0         1 

' 

'l         0 

0 

^ 

= 

0    cos  <^ 

-sin  * 

0     sine}) 

cos  (j) 

(7) 
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III.   IMAGE  RESOLUTION 

The  three-dimensional  object  can  be  considered  to  be  an  aircraft, 
with  maximum  orthogonal  dimensions  on  the  focal  plane  of  Xmax,  Ymax, 
Zmax.  These  dimensions  are  illustrated  on  the  top  and  side  views  of 
Fig.  2.  Since  the  image  is  to  be  digitaized,  the  discrete  maximum 
dimensions  can  be  given  in  terms  of  pixel  sampling  intervals,  x  and  y. 
A  pixel  is  defined  as  the  smallest  unit  sampling  area  (x,y)  of  resolution 
for  the  digitized  image. 

Xmax(quantized)  =  XN  (8) 

Ymax(quantized)  =  YN  (9) 

Zmax(quantized)  =  ZN  (10) 

These  dimensions  are  illustrated  in  Fig,  3.  The  quantization  of  the 
image  would  have  a  minimum  rotation  that  could  be  detected  by  a  1  pixel 
change  at  extreme  points  in  the  image  and  would  be  equivalent  to  the 
following  angle: 


tan  e  =  -j^  Cl ) 


This  relationship  is  illustrated  in  Fig.  4.  KN  would  be  equal  to  XN,  YN, 
or  ZN  depending  on  the  angle  of  rotation. 
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Figure  2 
Maximum  Imaqe  Dimensions,  Orthogonal  Views 
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□  =  1  pixel 


YN  pixels 


ZN  pixels 


Figure  3 
Discrete,  Orthogonal  Views 
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1  pixel 


Figure  4 
Minimum  Detectable  Rotation,  Orthogonal  View 


If  the  axis  of  rotation  was  not  orthogonal  to  the  focal  plane,  the 
minimum  detectable  rotation  of  an  equal  sized  object  would  be  much 
larger,  as  indicated  in  the  top  view  of  Fig,  5. 


cos 


min 


KN-1 


KN 


(12) 


KN     = 


1-cos 


(13) 


min 


Sample  values  of  KN  from  equations   (11)  and   (13)  for  minimum  detectable 
rotation  are  given  in  Table  1. 
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Figure  5.  Minimum  Detectable  Rotation 

Rotation  Axis  in  Viewing  Plane 
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TABLE  1 

Axis  Orthogonal  Axis  Parallel 
to  image  plane  (11)  to  image  plane  (13) 

(degrees)      KN  (pixels)  KN  (pixels) 

1            58  6566 

5            12  263 

10             6  56 


If  it  is  assumed  that  the  rotational  angles  are  measured  by  lines 
through  the  two  extreme  pixels,  the  geometry  and  error  sensitivity  would 
be  as  shown  in  Fig.  6.  For  a  maximum  length  of  KN,  and  a  1  pixel  rota- 
tion: 


1 

tan  6  =  (11) 


KN 


For  horizontal  errors  H: 


1 

=  arctan  r (14) 

KN  -  H 


1 

de  =  d  arctan  ^^ (15) 

KN  -  H 


For  a  change  of  variables 


1 

Til  =  (16) 

KN  ±H 
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dH 

dm  = ^  (17) 

(KN  t  H)2 


de  =  -^-  dm  (18) 

1  +  m^ 

dH 

de  =  ; 2—  (19) 

(KN  -  H)>1 


For  example,  if  KN  =  58,  H  =  1,  then  de/dH  =  0.0003077  radians/pixel  = 
0.01763  degrees/pixel.  This  sensitivity  for  horizontal  errors  is  very 
small ,  as  expected. 

For  vertical  errors: 

1  -  V  (20) 


6   =  arctan 

KN 


1  tv 

m  =  (21) 

KN 

dm  =  -^  (22) 

KN 


dV  KN 

de  =  5 T 5-  (23) 

KN^  +  (1  -  V)'^ 

For  example,  if  V  =  1  pixel,  and  KN  =  58  pixels,  de/dV  =  0.01722  radians/ 
pixel  =  0.99  degrees/pixel.  This  is  much  larger  than  the  horizontal 
error,  as  might  be  expected  by  observing  the  geometry  of  Fig.  6. 
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Figure  6. 
One-Dimensional  Rotation  Errors 
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This  section  has  derived  equations  (11)  and  (12)  to  predict  the 
minimum  detectable  rotation  of  a  three-dimensional  object.  These 
equations  apply  to  a  three-dimensional  object  with  a  two-dimensional 
maximum  projected  dimension  of  KN  pixels.  Equation  (11)  applies  if 
the  axis  of  rotation  of  the  three-dimensional  object  is  parallel 
with  the  camera  axis.  Equation  (13)  applies  if  the  axis  of  rotation 
of  the  three-dimensional  object  is  orthogonal  to  the  camera  axis. 
Equations  (19)  and  (23)  predict  the  rotational  errors  resulting  from 
horizontal  and  vertical  measurement  errors  in  extreme  points  of  the 
two-dimensional  image. 
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IV.   IMAGE  FILTERING 


Two-dimensional  matched  filtering  can  be  considered  to  be  an 
extension  of  one-dimensional  correlation  filtering  [5].  The  out- 
put of  the  two-dimensional  matched  filter  will  not  be  a  transformed 
image  with  some  type  of  improvement  or  enhancement,  but  a  correlation 
plane  whose  amplitude  at  a  given  location  corresponds  to  the  degree 
of  correlation  of  the  input  image  (signal  +  noise)  with  the  desired 
image  (signal).  This  is  equivalent  to  sliding  one  photographic 
transparency  across  another  and  adding  up  each  black  point  that  was 
aligned  with  another  black  point  on  the  second  image,  each  white 
point  that  was  aligned  with  another  white  point  on  the  second  image, 
and  each  grey  point  that  was  aligned  with  another  grey  point  of 
equal  shade  on  the  second  image. 

Several  practical  factors  must  be  considered  in  the  actual 
simulation  of  the  matched  filter.  The  input  image  and  the  refer- 
ence or  stored  image  could  both  be  compared  (correlated)  as  contin- 
uous grey  scale  images,  but  this  would  be  a  tremendous  storage  and 
computational  task.  Therefore,  the  reference  image  in  the  computer 
simulation  is  digitized  and  stored  as  discrete  X,Y,Z  points  (FACIT 
Model)  in  a  three-dimensional  coordinate  system  [6].  It  is  recon- 
structed, after  rotation  and  translation,  as  unit  amplitude  lines 
projected  on  the  focal  plane.  It  is  also  assumed  that  f  <<  y,  and 
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f  and  y  of  equations  (1)  and  (2)  are  known.  This  knowledge  of 
focal  length  and  target  range  yields  a  constant  image  scaling  fac- 
tor (k).  This  type  of  projection  is  illustrated  in  Fig.  7.  An 
example  of  an  aircraft  FACIT  Model  is  given  in  Fig.  8. 

The  input  grey-level  image  is  also  transformed  into  a  similar 
format  by  use  of  an  edge  detector  to  extract  outline  information 
[7].'  The  edge  of  the  image  is  defined  as  the  boundary  between  the 
two  regions  of  different  grey  level. 

The  reference  image  and  the  transformed  input  image  are  applied 
to  a  matched  filter  to  determine  the  degree  of  correlation.  An 
ideal  correlation  function  C(x,y,k,e,'j;,(|))  would  have  an  impulse 
response  that  could  be  modeled  by  delta,  (sin  x/x),  triangle  or 
exponential  functions.  The  impulse  would  occur  when  the  reference 
and  transformed  images  were: 

1.  aligned  on  the  image  plane  (x,y), 

2.  equally  scaled  (k)  and 

3.  at  the  same  attitude  {Q»^a)- 

However,  due  to  the  presence  of  noise  and  the  approximation  of  the 
three-dimensional  object  with  a  finite  element  model,  the  correla- 
tion response  is  not  delta,  (sin  x)/x,  triangle  or  exponential  in 
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Figure  7. 
Three-Dimension  to  Two-Dimension  Projection 
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Figure  8.  FACIT  Model  of  F-102 


24 


shape.   It  is  not  necessary  to  know  the  exact  correlation  response 
function  or  even  the  constant  component  due  to  a  nonzero  mean  value 
of  the  reference  and  input  image  variable  function.   It  is  only 
necessary  to  know  the  values  of  the  variables  x,  y,  k,  e  ,  ,|^  and  (j) 
at  peak  correlation. 

The  following  simplifying  assumptions  are  made  for  an  attitude 
{Q,i','t>)   only  correlation: 

1.  the  images  have  been  scaled  to  an  equal  size  (k) 

2.  the  peak  correlation  on  the  image  plane  will  indicate 
translational  alignment  (x,y) 

Fig.  9  illustrates  a  typical  x-y  correlation  between  a  noisy  image 
and  computer-generated  reference  image.  The  peak  correlation  at 
(0.0,  0.0)  pixels  indicates  that  both  images  are  aligned  in  the 
x-y  plane.  This  is  because  the  noisy  image  was  generated  from  a 
computer  reference  image. 

The  attitude  of  the  input  image  is  unknown.  Therefore,  the 
measurement  of  the  attitude  of  the  input  image  is  equivalent  to 
applying  the  transformed  image  to  a  bank  of  matched  filters  with  each 
"tuned"  to  a  specific  attitude.  This  process  is  illustrated  in 
Fig.  10. 
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X-Y    CORRELATION    PLANE 


Figure  9 

Sample  64  X  64  Pixel  X-Y  Correlation  Plane  Obtained 
from  32  X  32  Pixel  Image  and  32  X  32  Pixel  Reference 
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Figure  10.  Image  Processing  Flow  Diagram 


If  there  was  a  priori  knowledge  of  azimuth  e  ,  elevation  i>   ,  and 
roll  <p  y   to  within  10  degrees  each,  and  a  1  degree  resolution  of 
measurement  was  required,  a  brute  force  approach  would  require 
10  x  10  X  10  =  1000  correlations  or  1000  matched  filters.  Three  axis 
quadratic  interpolation  is  used  to  reduce  the  number  of  correlations 
required  for  the  given  resolution  of  measurement  [8]. 

It  is  assumed  that  the  variables  azimuth,  elevation,  and  roll 
are  uncoupled  in  the  joint  correlation  function  C(x,y,e,iK  4>).  Assum- 
ing that  there  is  a  near  orthogonal  relationship  allows  C(x,y,6  ,i|j  ,4) ) 
to  be  sequentially  maximized  (minimized)  with  the  independent  vari- 
ables: X,  y,  azimuth,  elevation,  and  roll. 
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The  pitch  correlation  function  is  modeled  as  a  quadratic  func- 
tion C(0): 


2 

C(0)  =  a  e  +  b  e  +  c  (24) 


3?   =  2  a  8  +  b  (25) 


Equation  (25)  is  set  to  zero  to  maximize  (minimize)  C(o) 


5    =     b  (26) 

max     —rr— —  ^^' 

c.   a 


Si  =  0  degrees,  62  =  8^  +  A6  degrees,  and  83  =  0i  -  Ae  degrees  are 
substituted  in  equation  (24)  to  solve  for  the  unknowns  a,  b,  and  c. 
The  solution  of  equation  (26)  is: 


=      (C(e,)  -  C(e.))  Ae  ,  . 

max        26(62)  +  20(83)  -  4C(ei)  ^'^' ' 


Similar  equations  are  written  for  C(ii>),  and  C(-4)). 

The  computer  simulations  listed  in  Appendix  B  implement  equa- 
ion  (27).  The  three  correlations  C(8i),  0(62)  and  0(83)  are  calculated 
with  Ae  =  5  degrees  for  a  rough  approximation  of  9  j,-  Then  three 
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correlations  C(ei  ),  C(92  )  and  0(63  )  are  calculated  with  as  =  2 
degrees  for  a  sharper  quadratic  approximation  of  9^^^^  .  In  addition, 
the  +2  and  j^5  degree  intervals  are  tested  to  make  sure  that  the  mid- 
point correlation  C(9i)  was  greater  than  or  equal  to  the  end  point 
correlations  0(82)  and  0(63).   If  either  is  not,  the  interval  is  shifted 
by  2  or  5  degrees  in  the  direction  of  the  larger  correlation  and  the 
mid-point  is  tested  again.  This  dual  pass  computation  reduces  the 
unwanted  coupling  effect  of  the  attitude  variables.  The  +^5  degree 
calculation  detects  the  general  angle  of  maximum  correlation.  The 
+2   degree  calculation  reduces  the  estimated  angle  error  by  not  being 
as  sensitive  to  unsymmetry  and  nonl inearities  in  the  correlation  func- 
tion. 

A  sample  of  the  dual  pass  azimuth  estimation  is  illustrated  in 
Fig.  11.  For  this  sample,  noisy  edge  data  of  a  simulated  aircraft  at 
an  azimuth  angle  of  -45.0  degrees  is  used  as  input  data.  This  data 
is  correlated  with  computer-generated  edge  data  from  a  stored  air- 
craft model.  The  initial  estimate  of  attitude  is  -48.0  degrees.  The 
final  azimuth  calculated  is  -44.11  degrees,  a  0.89  degree  error  from 
the  true  -45.0  degree  aircraft  azimuth  attitude.  This  is  a  typical 
maximum  error.  For  example,  if  the  initial  estimate  of  azimuth  is 
-42.0  degrees,  the  final  error  would  be  0.0  degrees. 
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This  section  has  introduced  the  use  of  two-dimensional  correla- 
tion to  determine  the  attitude  of  a  three-dimensional  object.  Equation 
(24)  is  used  as  a  second  order  polynomial  model  for  the  correlation 
function  of  each  attitude  variable.  This  approximation  reduces  the 
number  of  computations  necessary  to  determine  the  peak  correlation  for 
each  attitude  variable.  This  quadratic  approximation  was  demonstrated 
with  a  computer  simulation  that  tested  for  the  maximum  correlation  of 
an  image  as  a  function  of  attitude.  The  simulated  image  was  previously 
generated  by  adding  edge  noise  to  a  computer-generated  image  at  a  known 
attitude  of  -45.0  degrees.  By  starting  the  correlation  algorithm  with 
a  -3.0  degree  offset,  the  final  error  was  found  to  be  0.89  degrees. 
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V.  MATCHED  FILTER  EQUATION 

It  is  well  known  that  one-dimension  correlation  can  be  calculated 
by  a  simple  multiplication  of  two  functions  that  are  Fourier  trans- 
formed [9]. 

c(t)  =  J  X  (t)  h(t+T)  di  =  x(t)  *  h(-t)  (28) 

where  x  is  the  dummy  variable  of  integration  and  *  is  used  to 
symbolize  the  convolution  integral. 

C(w)  =  X(w)  H*(w)  =  X*(w)  H(w)  (29) 

where  w  =  Zirf  =  frequency  in  radians  and  *  is  used  to  indicate  com- 
plex conjugation;  that  is,  if 

H(w)  =  R(w)  +  jl(w),  then  H*(w)  =  R(w)  -  jl(w). 

c(t)   <=>  C(w)  (30) 

where  <=>  is  used  to  indicate  a  Fourier  transform  pair;  that  is 

C(w)  is  the  Fourier  transform  of  c(t). 
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This  same  relationship  can  be  shown  in  two-dimensions  and  in  discrete 
form  [5]. 

N-1  M-1 
c(l<,l)  =   Z    z   X  (i,j)  h(l<+i,l+j)  (31) 

i=0  j=0 

C(m,n)  =  X*(m,n)  H(m,n)  (32) 


K-1  L-1 

,•  2  nmk      ,•  2  nnl 
C(m,n)  =       z    z   x(k,l)e  '-'   K      e  -'  L 
k=0  1=0 


I-l  J-1 

.  2n  mi      .  2  n  nj 
z    z   h(i,j)e  ^       I  e  ^  ~^ 

1=0  j=0  (33) 


m  =  variable  in  first  dimension 


n  =  variable  in  second  dimension 


K  =  I  =  M  =  Period  in  first  dimension 


L  =  J  =  N  =  Period  in  second  dimension 
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A  complete  listing  of  the  transformation  indicated  in  equations  (32) 
and  (33)  is  given  in  Appendix  A. 

The  reason  the  Fourier  domain  correlation  approach  is  used  is 
that  for  large  data  arrays  the  Fast  Fourier  Transform  (FFT)  computa- 
tion time  is  less  than  that  of  direct  time  or  space  domain  correlation. 
Also,  Fourier  transform  and  array  processing  hardware  is  available  for 
many  scientific  computers. 

The  continuous  correlation  equation  (28)  states  that  an  infinite 
integration  is  necessary  to  determine  the  correlation  at  a  particular 
time.  This  is  not  the  case  for  finite  sampled  functions.  The  summa- 
tion interval  for  the  sampled  functions  requires  that  the  majority  of 
the  information  in  the  two  functions  be  contained  in  the  finite  sampled 
interval  for  the  correlation  to  approximate  the  continuous  correlation. 
When  x(t)  and  h(t)  are  sampled  with  an  impulse  function  A(t)  over  a 

finite  interval,  0<t<t    =  NT,  their  Fourier  transforms  x(w)*a(w) 

max  • 

and  H(w)*a(w)  will  be  periodic  with  period  w  =  1/T,  where  T  =  sampling 
interval  of  the  impulse  function  A(t).  When  the  continuous  function 
x(t)  is  sampled  over  the  interval  0<  t<t^^„  and  Fourier  transformed, 

niaX 

this  is  equivalent  to  calculating  the  Fourier  Series  coefficients  of  a 
sampled  periodic  sequence  x  (nT)  of  period  N,  where  N  =  t  /T  and 
T  =  time  interval  between  samples.   If  x  (nT)  is  time  shifted  in 
relation  to  x(nT)  by  m  samples,     and  the  first  N  samples  are  trans- 
formed, it  can  be  shown  that  the  resulting  complex  Fourier  transform 
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will  change  in  phase  but  not  in  amplitude  or  period  [9].  In  the 
time  domain,  this  shift  of  m  samples  is  due  to  the  periodicity 
of  X  (nT)  and  is  referred  to  as  a  circular  shift  of  a  sequence  [5]. 
The  sequence  appears  to  be  circular  because  the  last  sample  of  x(nT) 
that  is  shifted  beyond  the  interval  N  T  appears  at  the  beginning  of 
the  sequence.  Likewise,  in  the  two-dimensional  correlation  there  is  a 
circular  shift  in  two  dimensions  when  one  image  x(i,j)  is  correlated 
with  another  h(i,j).  The  two-dimensional  periodic  or  circular  sampled 
function  is  equivalent  to  joining  parallel  horizontal  and  vertical  edges 
of  each  image  and  rotating  one  image  on  each  of  its  axes  while  the 
other  image  remains  stationary.  This  two-dimensional  circular  shifting 
is  depicted  in  Fig.  12.   In  this  figure,  image  2  has  shifted  n  samples 
in  the  horizontal  direction  with  respect  to  the  fixed  image  1.   Image  2 
has  shifted  m  samples  in  the  vertical  direction  with  respect  to  the 
fixed  image  1.  To  clarify  that  image  2  has  circularly  "wrapped  around" 
image  1,  the  corners  are  labeled  A,  B,  C  and  D. 
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Figure  12.   Implementation  of  Two-Dimensional  Correlation 
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The  continuous  images  1  and  2  are  sampled  at  N  and  M  discrete 
points  in  the  x  and  y  axis  respectively.  This  sampling  results  in 
a  periodic  equation  (33)  in  the  transformed  frequency  domain  of  per- 
iod N  and  M.   If  the  correlation  product  of  the  two  sampled  images 
has  nonzero  terms  beyond  the  N  and  M  dimensions,  there  will  be  an 
overlap  of  nonzero  samples  in  the  frequency  domain.  This  overlap 
phenomenon  of  a  high  frequency  component  taking  on  the  identity  of 
a  lower  frequency  component  is  referred  to  as  frequency  aliasing. 
Likewise,  spatial  aliasing  would  occur  if  there  was  an  overlap  of  the 
periodic  spatial  functions;  for  example,  if  the  first  and  last  samples 
added. 

To  test  the  correlation  method  of  attitude  measurement,  a  com- 
puter simulation  program,  FAC  listed  in  Appendix  B,  was  run.  Discrete 
boundary  points  were  generated  as  input  data  for  a  known  aircraft 
attitude  as  illustrated  in  Fig.  13.  A  noisy  test  object  was  generated 
that  was  within  +_5  degrees  in  azimuth,  elevation  and  roll  of  the  ref- 
erence object.  The  noisy  boundary,  as  described  in  the  next  paragraph, 
was  compared  with  the  noise-free  computer-generated  edge  data  to  esti- 
mate the  aircraft  attitude.  The  representative  errors  of  attitude 
(azimuth,  elevation  and  roll)  were  computed  by  comparing  the  attitude 
that  was  iteratively  predicted  by  equation  (27)  with  the  actual  atti- 
tude. These  results  are  summarized  in  Table  II. 
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Figure  13 
Boundary  Points  of  FACIT  Model 
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Table  II. 


Representative  Attitude  Errors  (Degrees) 
Three  Sample  Average,  1  Pixel  Noise  Standard  Deviation 


Azimuth 

0.12 

Elevation 

0.55 

Roll 

0.71 

It  was  noted  that  the  errors  increased  with  each  attitude  measurement. 
This  is  because  of  the  sequential  coupling  of  a  previous  attitude 
angle  error  to  the  next  quadratic  angle  interpolation. 

The  quadratic  correlation  attitude  estimator  (27)  was  initial- 
ized with  data  to  within  10  degrees  of  the  known  aircraft  attitude. 
Boundary  points  on  a  32  x  32  grid  were  generated  from  the  two- 
dimensional  projection  of  the  test  aircraft.  The  number  of  boundary 
points  generated  were  a  function  of  the  test  aircraft  attitude,  but 
typically  108  points  were  generated: 

(27  each:  maximum  horizontal,  minimum  horizontal,  maximum 
vertical,  minimum  vertical). 


38 


To  simulate  detector  noise,  Gaussian  noise  with  zero  mean  and  one 
pixel  standard  deviation  was  added  to  each  boundary  point,  horizon- 
tally and  vertically.  Since  the  boundary  point  function  was  on  a 
discrete  32  x  32  grid,  this  had  the  effect  of  quantizing  the  noise 
to  unit  pixel  intervals.  For  example,  if  the  absolute  value  of  the 
Gaussian  noise  was  less  than  0.5  pixel,  no  noise  would  be  added.   If 
the  absolute  value  of  the  Gaussian  noise  was  between  0.5  and  1.5 
pixels,  one  pixel  of  noise  was  added.  Likewise,  two  or  more  pixels 
of  noise  would  be  added  for  greater  outputs  of  the  generator.  These 
noisy  boundary  points  are  shown  in  Fig.  14. 

An  uncalibrated  system  test  with  real  data  was  made  from  a 
digitized  video  image  of  a  F102  aircraft.  The  intensity  data  from 
the  image  was  transformed  into  edge  points.  The  edge  points  were 
correlated  with  the  model  data  stored  in  the  computer  memory.  This 
predicted  the  aircraft's  attitude. 

The  x,y  resolution  of  the  digitization  was  340  pixel  columns 
X  220  pixel  rows  and  the  intensity  resolution  was  1  to  256.  As  the 
photograph  in  Fig.  15  shows,  the  aircraft  image  makes  up  only  a  small 
portion  of  the  video  image.  Therefore,  a  64  x  64  pixel  window  shown 
in  Fig.  16  was  contrast  enhanced,  in  that  the  original  intensity 
range  of  103  to  156  was  linearly  transformed  to  a  1  to  256  intensity 
level.   It  is  noted  that  the  digitized  image  is  not  displayed  at  the 
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Figure  14.  Boundary  Points  of  FACIT  Model  with  Added  Noise 
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Figure  15 
Photograph  of  F102  Aircraft 
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Figure  16 
Digitized  64  X  64  Pixel  Image  of  F102  Aircraft 


42 


correct  aspect  ratio  due  to  the  unequal  sampling  ratios  taken  verti- 
cally and  horizontally. 

The  edge  data  from  Rosenfeld's  Detector  [7]  was  applied  to  a 
variable  threshold  iteration  that  was  adjusted  until  greater  than 
200  pixels  were  above  the  threshold.  This  adjustment  was  necessary 
so  that  low  contrast  edge  pixels  typically  at  wing,  nose  and  tail 
extremities  would  exceed  the  threshold.  These  extremities  are  most 
important  in  obtaining  a  given  attitude  accuracy.  Fig.  17  shows  the 
results  of  the  best  computer-generated  match  of  the  stored  image 
superimposed  on  the  edge  data  obtained  from  the  video  image.  The 
computer  predicted  the  attitude  to  be: 

Azimuth  =  -115.0  degrees 
Elevation  =  45.0  degrees 
Roll     =  -  5.0  degrees 

There  was  no  other  data  available  to  establish  the  attitude  of  the 
F102  aircraft  at  the  instant  the  camera  recorded  the  image.   It  is 
presumed  that  further  investigation  would  show  this  attitude  to  be 
correct  to  within  the  limits  of  the  image  resolution  and  the  edge 
noise  present. 
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Figure  17 
Best  Fit  of  Stored  Image  Pixels  to  Video  Outline  Pixels  of  Fig.   16 
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The  lack  of  sharpness  in  the  transformed  edge  data  shown  in 
Fig.  17  can  be  attributed  to  several  factors; 

1.  Low  x,y  resolution  (64  x  64). 

2.  Low  intensity  range  (103  -  156). 

3.  Uneven  lighting  and  shading  on  edges. 

4.  A  simple  edge  detector. 

The  last  factor  would  be  the  only  one  that  could  be  significantly 
changed  given  a  free-flying  aircraft,  fixed  camera  locations,  and 
natural  lighting.  A  review  of  methods  used  for  extracting  features 
based  on  edges  and  contours  is  included  in  a  survey  by  Lavine  [11]. 
More  recent  theories  and  experiments  on  edge  detectors  have  been 
demonstrated  by  Abdou  and  Pratt  [12]. 

This  section  has  introduced  the  use  of  two-dimensional  Fourier 
domain  correlation  as  described  by  equation  (32).  Representative 
results  of  three  computer  correlation  simulations  were  listed.  The 
simulation  compared  a  computer  generated  edge  model  to  a  previous 
computer  generated  model  with  added  edge  noise.  Finally,  a  sample 
correlation  of  the  computer  generated  edge  model  with  the  edge  data 
detected  from  a  real  video  image  of  an  aircraft  at  an  unknown  atti- 
tude was  described. 
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VI.  AIDED  AIRCRAFT  TRACKING 

An  example  of  the  utility  of  the  attitude  measurement  is  demon- 
strated by  utilizing  an  aircraft  positional  estimating  Kalman  Filter 
with  combined  radar  and  electro-optic  sensor  measurements  [10].  It 
is  difficult  for  the  radar-only  estimator  to  keep  a  close  target  in 
track  if  it  is  highly  maneuverable.  The  criterion  of  "highly  man- 
euverable"  for  a  particular  radar  is  dependent  on  the  following: 

1.  The  sampling  interval  of  the  radar. 

2.  The  range  gate  period,  or  the  sampling  aperture. 

3.  The  receiver  sensitivity. 

4.  The  return  signal  to  noise  ratio. 

5.  The  mount  azimuth  and  elevation  dynamics,  or  time  constants. 

6.  The  measurement  accuracy  of  azimuth,  elevation  and  range. 

7.  The  point  source  tracking  of  a  finite  length  target. 

These  parameters  and  others  contribute  to  the  ability  of  the  radar 
to  track  with  minimum  error.  By  modeling  the  aircraft  target  to  be 
driven  by  an  acceleration  function  that  is  dependent  upon  aircraft 
attitude,  more  accurate  range,  range  rate  and  range  acceleration  state 
estimates  are  possible.  The  acceleration  function  is  based  on  the 
well-known  aerodynamic  coupling  between  the  target  attitude  and  the 
direction  of  the  target  acceleration.  This  coupling  is  summarized 
in  the  following  comments: 
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1.  Major  accelerations  (lift)  are  normal  to  the  velocity  vector 

and  the  pitch  axis. 

2.  Positive  lift  is  more  likely  than  negative  lift  due  to  pilot 

physiological  factors  and  to  structural  loading  design. 

3.  Accelerations  in  the  velocity  direction  (drag/thrust)  are 

generally  smaller  in  magnitude  and  of  shorter  duration  than 
the  lift  accelerations. 

A  complete  multiple  measurement  system  using  the  previously  discussed 
imaging  attitude  measurement  concept  as  a  radar  tracking  aid  to  the 
positional  estimating  Kalman  Filter  is  shown  in  Fig.  18.  In  the  upper 
left  corner  of  Fig.  18  there  are  n  radars  that  provide  range,  azimuth 
and  elevation  measurements  to  the  Iterated  Extended  Kalman  target 
estimator.   In  the  lower  left  corner  there  are  m  video  cameras  that 
input  data  to  m  correlators.  These  correlators  first  scale  the  stored 
model  to  a  size  equal  to  the  predicted  input  image  size.  The  scaling 
is  a  function  of  target. range  and  system  gains  such  as  telescope  mag- 
nification and  image  pixel  sampling  area.  The  correlators  then  apply 
equation  (33)  as  described  in  the  previous  section  for  the  computer 
simulation  FAC.  The  attitude  measurements  azimuth,  elevation,  and 
roll,  from  the  correlators,  are  transformed  into  a  common  coordinate 
system.  The  predicted  attitude  of  a  model  aircraft  in  a  coordinated 
turn  is  combined  with  the  m  attitude  measurements  in  a  Kalman  Filter. 
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Figure  18.  Attitude  Measurement  and  Tracking  System 
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The  estimated  attitude  from  the  Kalman  Filter  predicts  an  accelera- 
tion vector  that  is  used  as  an  input  to  the  target  estimator. 

The  simulation  of  a  target  range  (R)  estimation  was  performed 
with  the  program  KAL  listed  in  Appendix  B.  The  simulation  KAL  was  a 
limited  example  of  the  radar  tracking  filter  as  illustrated  in  Fig.  18. 
This  simplification  was  made  possible  by  not  including  some  tracking 
variables,  such  as  azimuth  and  elevation,  to  the  full  three-dimensional 
trajectory  filter.  The  trajectory  was  truly  three-dimensional,  hence 
there  were  changes  in  azimuth  and  elevation,  but  they  were  not  used 
as  system  measurements  or  estimated  as  parameters. 

The  Iterated  Extended  Kalman  Filtering  subroutine  FTKALI  is  a 
modification  of  the  IMSL  Kalman  Filtering  subroutine  FKALM  [14]. 
The  system  was  modeled  as  a  radar  track  of  a  nominal  constant  accel- 
eration target.  The  assumptions  for  this  trajectory  are: 

1.  The  average  acceleration  rate  is  zero. 

2.  The  acceleration  rate  between  radar  scans  is  constant. 

3.  The  acceleration  between  scanning  intervals  in  uncorrected. 

The  above  assumptions  are  described  in  the  well-known  Alpha-Beta- 
Gamma  Filter  equations  [13].  Written  in  a  discrete  form,  these 
equations  are: 
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^  =  Vl  -^Tx^.^  +   2  Vl^  6 


n-1 


(34) 


\  -   Vl  ^  T  ^  n-1  ^  T-  ^"^n-l 


(35) 


•  •   • ' 


\   =  ^n-1  ^  "^  ^^n-1 


(36) 


where  ar  is  a  random  sequence  with  zero-mean  and  a  mean-square 
acceleration  rate  of 


^[''i   n] 


a  for  j  =  k 
0  otherwise 


(37) 


The  target  trajectory  was  simulated  in  three  parts: 

1.  A  target  constant  closing  velocity  of  100  or  175 
yards/second  at  the  same  altitude  as  the  radar. 

2.  A  10-g  target  turn  in  the  azimuth  plane  of  the  radar. 

3.  A  10-g  target  turn  "up"  in  the  geometric  plane  perpendicular 
to  the  azimuth  plane  and  the  ground  projection  of  the  radar 
range  vector. 
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The  final,  or  third  turn,  starts  when  the  target  has  rotated  through 
90.0  degrees  of  the  10-g  turn  in  the  azimuth  plane.  The  maneuver 
was  started  at  three  initial  ranges: 


1.  Distant 

2.  Mid-range 

3.  Close-in 


11275.0  yards 

1000.0  yards 

150.0  yards 


Inputs  to  the  filter  of  two  noisy  radar  range  measurements  and  an 
optional  radial  component  of  the  predicted  normal  acceleration  vector 
were  generated  at  a  sample  rate  of  10/second.  A  one  sigma  standard 
deviation  Gaussian  noise  of  1.0  yard  was  added  to  the  two  range 
measurements.  These  inputs  were  input  to  the  range  filter. 

Typical  results  of  range  error  for  one  simulation  of  all  man- 
euvers are  listed  in  Table  III.  The  addition  of  the  predicted  normal 
.acceleration  measurement  to  the  radar  range  measurements  reduced  the 
peak  and  RMS  errors  for  all  maneuvers.  The  reduction  in  error  was 
most  dramatic  for  the  close-in  maneuver  because  of  the  range  reversal 
for  the  particular  azimuth  and  elevation  geometries.  At  the  mid- 
range  and  distant  range  trajectories  there  was  no  range  reversal 
so  the  additional  acceleration  measurement  was  less  helpful  in  cor- 
recting the  range  and  range  rate  estimates  In  addition,  the  high 
radar  sample  rates  and  the  averaging  effect  of  using  two  radar  measure- 
ments reduced  the  difference  in  error  between  the  radar-only  estimate 
and  the  combined  estimate. 
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Table  III. 

Range  Error  (RMS  Average/Peak  -  in  Yards) 

Maneuver  Combined  Measuremert      Two  Radars  Only 

(Normal  Acceleration 
and  Two  Radars) 

1.  Distant  0.76/3.22  0.99/4.82 

2.  Mid-Range  0.89/4.25  1.06/5.53 

3.  Close-In  4.83/7.56  10.28/12.02 


Fig.  19  shows  the  aircraft  10-g  turn  trajectory  plotted  in  three 
dimensions  for  a  175  yards/second  air  speed.  As  compared  to  the  100 
yards/second  air  speed  trajectory  shown  in  Fig.  20,  the  distance 
traveled  down  range  is  much  greater,  and  the  turning  radius  was  much 
larger.  However,  within  the  5  second  interval  of  flight,  the  air- 
craft in  Fig.  19  did  not  climb  as  high  as  that  in  Fig.  20.  This  is 
because  the  third  maneuver  (turn  up)  did  not  start  until  approximately 
3  seconds  versus  2  seconds  for  the  100  yards/second  trajectory. 
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Figure  19 

Three  Dimensional  Aircraft  Trajectory, 
Velocity  =  175  yards/second 
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Figure  20 

Three  Dimensional  Trajectory, 
Velocity  =  100  yards/second 
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In  Fig.  21  the  radial  error  and  the  measurement  noise  in  yards 
is  plotted  versus  time  for  the  close-in  maneuver.  The  radial  error 
is  equal  to  the  true  radial  distance  minus  the  estimated  radial  dis- 
tance in  yards.  Gaussian  noise  with  zero  mean  and  a  standard  deviation 
proportional  to  the  square  root  of  the  estimated  radial  range  variance 
was  the  measurement  noise  added  to  the  true  range.  For  this  example 
the  noise  has  a  1  sigma  standard  deviation  of  2.0  yards.  The  wery 
large  peak  error  at  0.8  seconds  is  due  to  the  aircraft  flying  yery 
nearly  "overhead".  This  then  reverses  the  sign  of  the  relative 
velocity  vector. 

The  velocity  reversal  is  shown  in  Fig.  22  where  X(2)  is  the  Kalman 
estimator  velocity  state,  and  VR2D  is  the  first  difference  of  range 
divided  by  the  sampling  interval.  The  slow  time  constant  of  X(2)  is 
due  to  the  lack  of  a  velocity  measurement.  A  more  accurate  velocity 
estimate  is  shown  in  Fig.  23  for  the  distant  turn  maneuver.  There  is 
no  velocity  reversal  in  this  trajectory. 

Fig.  24  shows  the  range  acceleration  state  estimate,  X(3),  and 
the  measured  range  acceleration,  Y(3).  The  two  accelerations  track 
quite  closely  and  describe  two  phenomenon: 

1.  An  initialization. 
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Range  Error  and  Measurement  Noise 
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Figure  22 

Estimated  Range  Velocity  and  Differential   Range  Divided 
by  the  Sampling  Interval,  Close-in  Maneuver 
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Figure  23 

Estimated  Range  Velicity  and  Differential  Range 
Divided  by  the  Sampling  Interval,  Distant  Maneuver 
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Figure  24 
Range  Acceleration  Estimate  and  Range  Acceleration  Measurement 
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2.  Cyclic  acceleration  matching  the  period  of  the  third 
part  of  the  maneuver. 

The  acceleration  state  estimate  in  Fig.  25  with  no  acceleration 
measurement  does  not  show  this  periodic  acceleration,  indicating  the 
uncertainty  present  in  the  state  estimate. 

This  section  applied  the  correlation  attitude  measurement  tech- 
nique to  the  radar  tracking  problem.  The  attitude  is  related  to 
acceleration  and  therefore  a  useful  measurement  in  predicting  target 
acceleration.  Recursive  trajectory  equations  were  used  in  a  computer 
simulation  that  estimated  aircraft  range.  The  simulation  demonstrated 
that  the  Kalman  Filter  was  more  accurate  in  tracking  a  target  with 
varying  trajectories  when  the  predicted  acceleration  was  used  as  in 
addition  to  the  noisy  range  measurements. 
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Figure  25 
Range  Acceleration  Estimate  Without  Range  Acceleration  Measurement 
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VII.   CONCLUSIONS 


Manual  optical  image  comparisons  have  been  the  primary  method 
of  measuring  attitude  of  three-dimensional  objects.  Section  V  lists 
the  general  two-dimensional  Fourier  transform  equations  that  are 
applied  to  an  images  edge  data  and  a  stored  model  of  the  solid  ob- 
ject. The  attitude  of  a  three-dimensional  object  is  estimated  by 
using  a  three-axis  quadratic  interpolation  to  sequentially  maximize 
the  correlation  output  as  a  function  of  the  attitude,  position  and 
size  variables.  Moreover,  this  measurement  or  estimate  is  equivalent 
to  the  maximum  output  of  a  bank  of  matched  filters. 

This  concept  was  then  applied  to  realistic  Kalman  Filter  track- 
ing simulations  and  shown  to  give  improvement  over  conventional  target 
state  estimators  that  do  not  use  the  attitude  measurement  as  a  pre- 
dictor of  target  acceleration. 

Since  the  necessary  attitude  measurement  is  within  the  electro- 
optical  tracker  and  computer  state-of-the-art,  the  use  of  this  esti- 
mator concept  can  improve  present  tracker  system  capabilities. 
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APPENDIX  A 
TWO-DIMENSIONAL  DISCRETE  CORRELATION 

This  appendix  demonstrates  that  the  two-dimensional  discrete 
correlation  transform  (31),  (32)  is  equivalent  in  both  the  time  domain 
and  the  frequency  domain. 

N-1  N-1 

c(k,l)  =   z   E    x(i,j)  h(k  +  1,  1  +j)        (31) 
i=0  j=0 


C(m,n)  =  X*(m,n)  H(m,n)  (32) 

where  *  =  complex  conjugate 
Equations  (31)  and  (32)  are  verified  by  performing  the  Fourier  Transform 
and  conjugation  operations  indicated. 


H-1  N-1 
c(k,l)  =   2   z 
i=0  j=0 


1 


.  2n  .    .  2n 
-J  -^1m   -J  -^  J  n 

X*(m,n)  e        e 


N 


M  N  m=0  n=0 


P-1  Q-1 

1  ■  z   z 


2n  2n 

j  p(k+i)    j  q  (1+j) 

H(p,q)  e     M  N 


PQ  p=0  q=0 


c(k,l)  = 


(A-1) 


MN 


M-1  N-1  P-1  Q-1 

E    Z 

m=0  n=0  p=0  q=0 


2n 


2n 


211 


j  -M— P  k  j  -^-^   1 
X*(m,n)  H(p,q)  e   "     e   " 


2n 


2n 


(A-2) 


2n 


N-1  N-1    -j  -jjj— i  m   j  -n— j  n   j  ^;p-p  i   j  ^-q  J 
1    z   z   e         e        e        e 


L  MN   i=0  j=0 
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The  term  in  brackets  in  equation  (A-2)  is  equal  to; 


MN 


(M  N)  =  1 ,  if  m  -  p  and  n  =  q 


0       ,  if  m  ?^  p  and  n  ?*  q 
due  to  orthogonality  of  i  and  j.  Therefore: 


1   M-1  N-1  ^^       .  ^^ 

c(k,l)  =-j;}|^   z  I     X*(m,n)  H(m.n)  e  ^'     "^  ^  -  ^'     "^ 
m=0  n=0 


M 


N 


=  2 


-1 


X*(m,n)  H(m,n) 


(A-3) 
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APPENDIX  B 
COMPUTER  PROGRAMS 

This  Appendix  is  composed  of  the  listings  of  the  Fortran  programs 
ind  subroutines  for  the  attitude  measurement  and  target  tracking 
emulations.  The  programs  were  run  on  a  Univac  1110  computer  with 
I  Fortran  compiler  and  a  FLECS  translator  (Fortran  Language  with 
Ixtended  Control  Structures)^  The  graphical  outputs  were  plotted 
ising  DISSPLA  programs  that  generated  COM-80  crt  images. 
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ill 

131 


15 


Q^  nf-. 

u  I  >■  c  N  b  i  f  N  W  H  (  1  c  2  0  ^  )  .  *  h  I,  (  1  0  3  .'  *.  ) 

U  I  •)  k  IM  i  I  0  N  ii  (  1  2  O  )  .  1  u  V  (  1  c     )  .  M  (  3  ) 

Ul.iLNSItN  AZA'&^,cLA(6>,hULLA(t^,Cr'AX^c") 

U  I  .^t  t  u  i,  I  0  N  Z  (  c:  "i  o  ) 

CO;V)h'LtX*OA(165oA) 

coi»l^o^    tehrt.wH.A.Cfi^x.z 

DO     111     J=1,2u 

C  1^  A  X  (  J  )  =  - 1  ,  0 1  1  0 

CON  1  INUc 

I  PL  cr  =  i 

I  P  K  1  N  T  =  1 

KFAb      Ci.1v.)      AZ.tL.ixOLL.LAbT 

F  0  i^  1-  A  T  (  ) 

».RITt(o,1j1)     AZ.cL. ROLL. LAST 

FOK!HT(1X,3Fy.2fI'.> 

IFCLmST.LQ.I)        Ih'LuT  =  :5 


READ     ARRAY     SIZE     IN     PO*EKS     OF     ^ 
HEAL     ^U^eER     OF     AWhAVS     TO    Bt»'.    COkRELATED 


READ(b,5)K<,jJ 
F  0  R  ^  A  T  <  ^  1  5  ) 
REA  D  (  i  .  5)     NA 
Nl  =  Z  *  "KK 

n3=  1 

KK2  =KK*  1 

J J2  =J J*1 

M2=£*N1 

N22=i:*^2 

NH  =  M/2 

NJ=N2/2* .5 

N  V  =  r.  1  *  N  i 

NV4=r^1<f*N<:2 

NC=  (M*Ni:)/  2*»M 

!«  ( 1  )=KK  2 

h"2  >  =  JJ  2 

N  (  3  )  =  0 

FI  =  5.U159 

WRI Tt (6.1 5)     Nl ,n2 

FOhhMK'  I.^PUT 


^MA  X    OF     20> 


ARRA^'     Dl^'tMblONb     =     (  '  ,  I  J^.'X' .  I3  ,  ">  '  f  /  ) 


HEAD      IN     THE      UMlA     AhRAY 

*  R  I  T  t  f  0  .  ii  > 

DO     n    k=1 .nc 

K  1=  (K -1  ) *  ^4+  1 

k24=n"2'« 

f<EAC<5,7)'*H(I),I  =  Kl,Kt'i) 

«^'I^b(6.'/)<1.(•*r(I).I=^'?.^2'i) 

CON  T  i  N  U  L 

F  0  H  r«  A  T  (  c  *,  F  ,'  .  '.  ) 
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9        F0t<l-AT(1X.I5.24(1X.F4.0)) 

.      OOUllE     Zllt     OF     DAT„     AKKAY 

wH     =     INPUT     \JtCTQh,      wHK     =    OUTPUT     VtCTOR.     NV     =     li^PUT     DINLNblO** 

CALL     VtiOUbL  ( -H  ,  »(HR  .  NV  .NV<.) 

PLO  T     DATA     ARkAT 

TRANSFER     THE     WttU)     DATA     TO     THE     COMPLEX     A     Mi^TKlX 

DO  i  12  =  1  .^2<i 
00  3  11=1  ^f4^^ 
IWH=(IZ-1)*Nl'»ir 
A (I WH)=WH« ( ImH) 
CONTINUE 

CO.IPUTE  THE  UlbCRETE  FOURIER  THANbFORM  OF   THE  IkH(K)  DATA 

CALL  HARM (A ,M,  INV  .i  ,1  ,  IPtRR  ) 

READ  IN  THE  SECOND  DATA  ARRAY 

ROTATE  ►'OOEL  IN  AZMUTh 

NL=C 
NA1  =1 
DELAL  =5. 
AZA  ( 1  )  =  AZ 
AZA (i)=AZ-DELAL 
AZA(i)  =  AZ-»     DELAL 
290  CONTINUE 

00  3CG  J=NA1 ,NA 

GENERATE  I»-AGE  dOUNQAKV  ^HQt*     3 -D  I  M  E  N  Tl  ON  AL  MOuEL 
wHR  =  CmRD  =  OUTPUT  VtCTOR  (INCYM  X  INCYM) 

CALL  SURFAC(A2A(J),EL.R0LL.YMAX,Zf«AX,Y»lIN.^.^IN.J.IPL0T.IPRI^T) 

CALCULATE     CORRELATION     BETWEEN     I»JPUT     IWAGE     AND     GENERATED     IMAbE     30UNDA 

2G37     FOR  NaT<  /  ,  5X  ,  '  J     =' ,  1  2,  :.x  ,  'CM  AX     =',F17.G,/) 

CALL     CALCOR  (KK2.  J  J  ..  .NV     .  AZ  A  (  J  )  .  E  L  .  K  OL  L  .  J  .  N  1  2  ♦  N2  2  ) 
3  0u     CONTINUE 

NL=NL ♦ 1 

IF(NL.GT.6)     \jO     TO     ^202 

WRITc(6,310) 
ilL     F0RI«'AT(/,5X.'J'.3X,'AZA(J)'.7X,'C^AX(J)  './  J 

DO     311     J=1.3 

wRITEU,312>     J,AZaU),CI^AX(J) 
3  11     CONTINUE 


1i  ft''^'''^T(HX.I2,3x,Fo,2,3x.Fn.i;) 

TtST  fir,     f^AXiMU.M  CORklLATION.   U   C1'Ax(1)  <   CrtAX(?)   OH  C^Ax^j) 
CuFkELATb  AT  ADUITIO.M^L  H  U  1 1.  T  S 

IF(Li'iAX(1).LT.CiViAa<-))  Gu  TC  il4l 
I  F  (  C  I"  A  X  (1  ^  .  L  T  .  C  /•  A  X  (  M  ^  (,  (J  TO  3  1  A  il 
bO  TO  i-ii,^ 
iUI  IF(CMAX(Z).LT.CMAXi3))  Gu  TO  iU^ 
aza  ^  3)  =  A2  aM  ) 
AZ A  (  1 )  =  A2  A{2  ) 
AZA  («:)  =  A2  A(  2)  -DfcLAL 
C>AX(3)=C^AX(1) 
Cf>^AX(1)=C^>AX(2) 
CfA  )i  («:)  =-1  .E  10 
NA  1  =c 
NA=;. 
GO  TO  2vO 

}'hj.   continue 

AZA  (2)  =  AZ AC  1  ) 
AZA ( 1 )=AZ A(3) 
AZA(3)  =  AZA(3)*Dt*LAL 
Cf«A  X  (2)  =L^  AX  (  1  ) 
C  w  A  X  (  1)  =  C  !^  A  X  C  3  ) 
C^»A  X  (3)  =-1.E10 
NAl  =i 
NA=  3 

0  0  To  c  V  G 
iU*.  CONTINUE 


CALCULATE     >)OST     PRUaABLt     AZMUTH     FOft     INPJT     FkAME     BY     «UDRaTIC     InTc^F 

AZST=(CHAX(2)-Li«<AXl3))*DELAL     /(2.*C.'«AX(I)*  i:,*CMAX(2)-4.*CMAX  (1  )) 
lAZbT=IFIx(AZST*.5)*IFIxrAZAri)> 
222c     F0R^'AT(1A6) 

PRINT     320,AZST.IAZ:>T 
320     F0.<^AT(/,5X,'AZST     =',Fl2.5,5x,'lAiST     ='.I4,/) 

KEDUCE      INTERVAL     OF     INTERPOLATION     if     AdS.     ^AlUc    OF     AZST     >     1. 

I  F  (  At:S(  AZ  ST)  .L  T  .1  .o)     bU     TO     3202 
DEL  AL  =  2  .0 
NAI  =2 
N  A=  3 

AZA  (2)  =  AZ  A(  1  )-uELA,. 
AZA  f3)  =  AZAn  )*DELAl 
CMA  X  (2)  =-1 .E 10 
C^AX (3) =-1 .E 10 
GO     TO     29u 
3^"  4     CONTINUE 

AZ=  FLOAT(  lAZST) 

KuTATt     i'iGuEL     IN     ElcVAIION 


<;p 


DELi.Hb  =  5  , 

N  L  =  :. 

C  :■*  A  X  (  7  ^  =  C  1*^  A  X  (1  ) 

NAiJ  =<i 

N  A«;i=v 

ELA  (l)  =  tL 
tLA<^)  =  EL-OELfc>o 
cLA (3)=tL*DELEPS 
:i9u     CONTINUE 

DO     UoC     J=NAo,NAV 
JK= J-C 

GtNtPATE     I^•AbE     uOUNOArtY     FkOM     3 -D  I  ri  E  N  T  I  ON  A  l     I^Oi/EL 
Whfi     =     CmKH     =    OUTPUT     VtCTOr<     (INCYM     X     INCY^) 


CALL     SUKFAC(A2,tLA(JM),ROLL.yMAX,Zf*AX,Yf^IN.2l«»lN.J.lPLOT.IPRINT) 
CALCULATE     COf<«£LATION     bETwEEN     IMPUT     IMASc     AND     GENfcRATcU     I M  A  vj  E     30UNDA 

CALL     CALCOR  (KK<:  .J  J2  ,N\/     .  AZ  .  EL  A  (  J  H  )  .  ROLL  .  J  .  M  2  .  ?^  2  ^  ) 
^GG     CONTINUE 

NL=NL ♦ 1 

IFruL.GT.6)     GO     TO     nidi 

WRI Tt (6.410) 
<,1'^     F0R^■AT(/,bX,'J'.3A.'£LA(J^«)',7X,'CfflAx(J)'./) 

00     411     J=7,9 

bRITt:(0.312)     J,tLACJM).CWAX(J) 
41  1     CON  TiNUE 

TEST    FOR     MAXIMUM     CORRtLATION.     IF     CMAXC?)     <     CHaX(8)     OR     CMAaCv) 
CORRELATE     AT     ADDITIONAL     POINTS 

I F( CMAX <7)  .L T  .CNAX (8) >     GO    TO     4141 
IF( CMAX (7).LT.CMAX(9) )     GO    TO     4142 
GO     TO    4144 
4141     IF  (  CtiAX  (6)  .L  T.Cl^AX  (9)  )     GO     TO     4142 
ELA  (3)  =  ELA( 1  ) 
ELA  (1)  =  ELA(2) 
ELA  (2)  =  EL  A(2  )-OELbt^S 
Cf^A  X  (9)  =CHAX  (  7) 
Cr*A  X  (7)  =C!*AX  Cb) 

cr^AX  (b)=-1  .ElO 

N  Ajj  =0 
NA9=0 
60     TO    390 
41c<:     CON  TINUE 

ELA  (2)  =  EL  A( 1 ) 

ELA  (  1  )  =  EL A(3) 

ELA  (3)  =  tLA(3 ) ^uELEPS 

C»"A  X  (b)  =Cf^AX  (  7) 

Cf*A  X{7)  =CMAX  (9) 

C^AX(9)=-1.E10 

N  A  V  =  9 


b  0     T  u     3  y  0 
-.l^H     CCNTiNUc 


CALCULATE     MOST    PROcABLt     cLEVATlON     FOR     INKUT     FPA^t     tY     UuDhATIC     I'il 

ELST  =  (CNAX(d)-CI-AX(9))*PtL6PS/(2.*Cr'AX(9)>i:.*CMA*fc>-A.*CMAX<1)) 
IELST  =  IFiX(ELST*.5>*IfIMELA(1)) 
PRINT     3i:1,£Li,T.IELwT 
221     F0R^AT(/.5x,'ELbT     ^  '  ,  f  1  2  .  5  .  Sx  .  '  I  E  LS  T     ='.1*../) 

KtDuCt      INTERVAL     OF     I N T L R F OL AT  I  ON     IF     AUS,     wALUt    OF     ELST     >     1. 

I F (  AdS( cLST)  ,L T  .1  ,u)     bO     TO     420^ 
DEL  fcPS  =  2.. 
NAli=d  , 

NA9  =  S' 

ELA  (2)  =  EL A(1  )-bELtPS 
ELA  (i)  =  ELA(  1  )-»UfcLtPS 
CN'AXC     d)  =  -1.ElC 
Cf^'A  A  (     9)  =  -1  ,t10 
bO     TO     390 
'.c^'^tCONTINUE 

E L = F L 0  A T (I E L S T ) 

ROTATE     WODEL     IN     ROLL 

DEL  RH0  =  5. 
NL=L 

C«AX(  13)  =  CMAX (  7> 
fiOL  LAC  1  )  =  ROLL 
ROLLA(2)=R0LL-0ELKhO 
ROLLA(3)=ROLL*DcLRnO 
NA1  A=U 
NA15=15 
49C    CONTINUE 

00  5C0  J  =  NA14  ,NA15 

GENERATE     IhAtfc     dOUNDAWY     FkOM     3 -0  I W E N TI ON Al     hOuEL 
yHR     =     CARD     =     OUTPUT     VtCTOH     (I  N  C  Y  ^<     X     I N  C  Y  l«  ) 

CALL     SUnF«C<AZ,EL.KOLLm(JM)  ,Y!«AX,ZMAX.YMIN  .ZMIN.J.IPLOT.IPkINT) 

CALCULATE     CURRELATION     BETWEEN     INPuT     IhAGE     AND     GENERATED     I«AoE     9OUN0A 

CALL     CALCCr?(M<2.JJi,NV     ,AZ.EL.KOLLA(Jr^).J,i>12.N22) 
500     CONTINUE 

NL  =  NL*  1 

IF(NL.bT.fc)     bO     TO     :3  20  2 

wRI Tt (6.3  10) 
510     FOR»*AT(/,5x.'j'.3x.'ROLLA(Jf)'.7X.'CMAX(J)  './) 

DO     511     J  =  13,  15 
J!*=  J-12 


t'f'ITtCo.Si:)      J.rOLLA(JM),CWAX(j) 

i>  1  1     C  C  N  T  1  r.  u  E 

TlST     fCK     MAXIMUM     CUKHlLATIOn.      IF     CCAX(13)     <     CrAX(14) 

COM^ELMTt     AT     AOoiriUfv.ML     HOli'iTS 


UK     C  h  A  X  (  1  5  ) 


5U1 
5U£ 


U1 


GO     TO     5  1  ^, 


I  F  (  CNAX  (1  i)  .LT  .Cr-iAA  di,  )  )     GO     TO 

IF((.MAX(1^).LT.Cr'.AA(15))     GO     TO 

bO     TO     iU^ 

I  F  (  C  fi  A  X  n  4  )  .  L  T  .  C  ^  A  A  ( 1  5  >  ^ 

RCLLA(i)  =  fi3LLA(1) 

KOLLA(l)=ROLLA(i) 

ROLLMf2>  =  ^0LLA(c)-UELf<hO 

CrMXn5)  =  CVAX(13) 

C  ^  A  X  (  1  3  )  =  C  «  A  X  (  1  ^  ) 

C^A  X  (  1<,  )  =  -1  .  t  lU        i 

N  A  1  A  =  U 

N  A  1  b  =  U 

GO  Tu  i,^ij 
-!<*<:    CONTlNut 

f»0LLA^2)  =  S0LLA('i) 

KOLLA (1 )=KOLL A ( 3) 

f<0LLA(5)  =  R0LLA(i)*uELHH0 

Cl»«AX(U)  =  ClAX(1i) 

CWA  X  (1  3  ^  =  C«A  A  Mi) 

CfMX(U)  =  -1.c1U 

NAI  A  =  1i) 

NA1 5=15 

GO  TO  490 
514^     CONTINUE 

CALCULATE     >'OST     PROaAtJLfc  ROLL     FOR     INPUT     FhAMt     S»     -UDRaTIC     iNTt^P 

ROLLST=(CfAX  (U)-L,^AX  (1S))*DtLRH0/(2.*C^'AX  l15)*2.*Cf^AX(1^) 

*-4,*c^'Ax(1)) 

IR0LST=IFIX(RuLLST*.S)*IFIX(R0LLA(     D) 
HRINT     323.ROLLi»T.     iROLST 
323     FOR^•AT(  /,  5X. 'RULL:>  r     =',F12,5,5X.'     IROLST     =',I4./) 

RECOCE      INTERVAL     OF      I N  1  E  R  P  OL  AT  I  Oi>J     IF     Aa  S  .     wALUc    OF     ROLLST     >     1. 

IF  (  AbS  (ROLLS  T  )  .LT  .  I  .0)     GO     TO     bil^Z 
0ELRh'O  =  t'. 
NAI  A=  H 
NAl 5=15 

R0LLA(2)=K0LLA(l)-ufcLHH0 
RCLLA(3)  =  K3LLA(  l)*^ELr»-0 
CMA  >  M  ^  )  =  -1  .  c  1U 
CMA  X (  15  )  =  -1  .Elu 

0  C  To  <,  V (; 
52 Cc     CONTINUE 

IPL  01  =  2 

1  F  (  L  AST  .t  ^  .^  )      GO     Tvj      1 
STOH 
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S  Lk  M  C 


bUwF  «C     ROTATES     A     SlLHED      i -0  I  N  c  i^  T  I  0  N  A  L     MODtL     I  N      5     A  x  I  i        A  2  >  U  T  m  .  h  L  E  V  A  T- 
iCi\      Tht     HKOJECTtU     ll"Ab6     OF     The     i'«ODcL     ON    A     f-LANt     SUnFnCt     IS     OUANTIZtO 
uMO     M     (32>3i:)     LKIO.       OUTPuT    A^kAY     =    'CARu' 


J  buh  OU  T  I  (iE 

IPhINT) 
P  Ah/^l'-t  Tb  R 
P  AriAh'fcTkR 
?  AKANE  Tt  R 
P  AwAf^E  Tt  R 
D   INt  NS  ION 


SbKFAC(AZ.cL,^^0LL,Tf■AX1.Z''AXl,]n^lN1.Z^iNl.JA.IPL0T, 


NC  S  =  71  ,K  1  =  2  "NCb 

I  N  C  Y  M        =  5  3 

IZ  E=C.  I£=1 

NI  =5  .No  =  0.N0P=7 

AK^AX(I,\CY^').P(3,n.NCS).T(3.5.NCS).XN0Rrt(uCS),YK<2.NCS). 
AY  ^0K^  (NC  S)  .Z  NORf*(NCS)  ♦  AKHIN  UNC  Y[*)  .  bKMAXdSCYM        ). 

t)      bKMlNUNCVM        ) 
D   If-cNi  ION     CA  RO  (163csA  ) 
0  IMENilON     YL  (  1Cj)  .ZL(  10) 
C  CMi^ON     CAPO 
C  cr«i«ON/  F  ACOA  T/     P 
2  ERU  =  »J.O 
0  N£=  1  . 
L  CTR=0 
D  RC  =  ACOb (-^ .  )/ loC  . 


L-ENhKAL     PREPERATIOS     OF     DATA 

A  2R= AZ  *LRC 
E  LR  =  EL*URC 
RC'LLH  =  RuLL*DRC 
S  AZ=S IN ( AZR) 
C  AZ  =  CUS (AZR) 
3  El=  SIN ( ELR) 
C  cL=COS (ELR) 
S  kO=SIN ( POLLR) 
C  hO=COS  CROLLR  ) 

Y  NAX=-1 U.EIO 

Y  M  N  =  1  Cj  .  E  1  n 

Z  NAX=-1U.E10 

z  r^  I  r.  =  1 0 .  F 1  r 
1^22    F  CR^AT(  1  X.  15  .  UF  10  .2) 

hUTAT  E    TARGE  T 

DGlO»<  =  1.NCb 

DO    1  G    J  =  1  .  A 

T  (1.J.K)=P(1  .J.K)*C£L*CAZ-PU^J.K)«CFL*SAZ*P(3.J.iC)*S£L 

T(2.J.K)=PM.J.K^*<CRu*iAZ*oRO-Sf;L»CAZ)*P^2.Ji<)*<C^O«CAZ-ShO*ScL* 

IS  AZ  )  -P(  ^  .J  .Ik  )  *SKG*CEL 
T(i.J.K)=P(1.J.N)«(SRO*SAZ-cK0*i.!rL*CAZ)*P(c.J.K)-(CnO»SEL*S«Z 

1         iS RO*C AZ ) ♦ P( 3 .  J .K )  •CRo  •Ctu 
I    F  (  r  r  2  ,  J  ,  K  )  .  o  7  .  V  ••  A  X  >      Y  M  A  X  =  T  (  <i  ,  J  ,  K  ) 
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T  Kr(..J.K>).LT.VMN)     yMIN  =  T(2.J.K) 


ftVfc      FIRST     YW--        AND 
L  CTH=LCTR*1 
I  F  (  J  A  .  fM  t  .  1  )     Of     T  C 

V  c  A  X 1  =  y  t"  A  X 

y  M  N  1  =  y  n  I  N 
7  p-HA  1  =  ZrAX 
Z  i-  I  N  1  =  Z  i>  I  N 
C  OM  luUc 


if^PU  Tc     SURFACE     NORMALS 


Z^i--     V  ALUc  S 
1£ 


!:L^^CoI;T<Ll...-T<e.2.K...<T(..1.K>-T<3...K»-CT.2.1.K,. 

1  T(3.A.»(  ))*(T(1,1  .K)-T(  1  .2.K)) 

ZNOHMK)  =  <T(1,1,.)-T(1.2,K))*(T(2,1.K)-T^2,..K))-(TC1,1.K) 

1  T(1.'..lC))*(T(2,1.K)-Tlw.2.K)) 

15C0NTINUt 

:0l*PU  TE     lYMAX     ANU     IZHAX 

L  CTk  =  LCTR*1 

R  ZKAX  =  AbS(Zf1AX-ZMN) 
R  yKAX  =  AbS(yMAX-yHlN> 
I  ZMAX  =  RZf*AX 

I  Y^.AX  =  Ry^*AX 


IHREPA  RE     PLOT     PA&E 


TO     19 


K 


1   F  (I  PLOT  .EQ  .  C)     oC 

C  ALL     TMHPLT( V)  ,  V        vvTK      W 

AXL.TH=     AK.AXKCZi'AX     -  Z  r' i  N      )  ,     ^  ^' AX     -YYIN      ^ 

AXLbTH     =     AlNT(AXLbTH/2.)     *      lU 

X  C:>.TR  =  AInT((  y^.AX♦y^'IN)/2) 

y  cnTW  =  A1NT((  iMAX»ZK'IN»>/2) 

X  Ck^3  =  xCNTR-A  XLbTh 

Y  Lfib  =  YCNTR-A  XLbTh 

X  tND  =xCr,TR*A  XLbTF 

y   tNt>=yCNTR*AXLbTH 

CALL     abUPL'-D 

C  ALL     NOd  ROR 

C  ALL     »-ElGHT(  .1) 

£NCUDE(ca.2lL.LML-EL)AZ,EL,RuLL 

CALL     TITLEaA.EL,lL^,'         .     .        /J) 

C  ^LL     i-lAhKER(  '   J 
C  ALL     SCLPIC  (  ...c) 
f  UM^  A  T  ^      ' 


AZ  =  '  ,Fc.l 


.'FL  =  '»F^-''«-^''^^'-'-"''^**'' 


'  > 


7.-^ 


TAhT      Ii^.L=K      FUR      z-KAoTtk 

«  i,C  Y  I'  =f  lOAT(  INC  Yr'  ) 

0  L  r  =  K  Y  !•■'  M  X  /  R  sj  C  Y  N 

u  L;A  =  HZl^Ax/f<^cy^ 

1        FcKf''AT('  PLT1.        PLTE, 'c:6X  .  'YMN.' ,  1CX  .  'ZrIN' .1  ^X  .  'DL  Y  .  '  .  IcX  . 

C       'uLZA',      ttX.'lZ^AA,',      1A.'IY^■AX',/) 
Z        F  CK>'AT(0F16.  5,  il  L  .  /) 

K    l»A>  =U 

K  p  1  r  J  =  u 
i  ecu  =1  . 

1  SEtU=1 

«  hlTfc  (6  ,  1ii22  )        LCTk 
D  C    1  Cl)     I  Z  =  1  .  I^CY^ 
Z  tAK  =  IZ«DLZA*ZMIN 


UNi)      Y    INT6KSECT  IONS 


AKMAX(IZ),Dy(«lN.ZYMIN(lNCY),DYMAX. 
IFLG  .  lAKf Lb  (iNC Y)  .  I N C Y ' . / ) 


L  Vt=C 

A  KKA  X(  I Z  >  =  -1  O.EIO 
A  Kf-IK(U  )  =  1&  .610 
251        FCS?«AT('        Y,     ZBAR,     AKf«IW(IZ). 
.Z   YTcNPfHAX(INCY).',2  0X,'lZ.i, 
DC    ««G     1  =  1, SCS 
Y  K(1  . I )  =  YMIN -lUO  . 
I   FCXNORhd  )  .LE  .0.  )     GO     TO     AJ 
N  fc  k  =  0 

DO    db     H  =  1  .  A 
K  E  =  K-»1 

I   F(R  .EQ  ,4)     K  6  =  1 

I   F(AbS<  T^3.K  ,I)-T(3,KE,I  >).lE.1  .E-6)     GO     TO    c5 
Y=(T(2.KE.I)-T(2.K.I))*CZeAK-T(3.K,I))/(T(i.Kc,n-T(3,K,I))* 

1  T(i:  .K  .1  ) 

X  =(T(1,KE,I)-T(1.K.I))*(Z6An-T(3.K,I))y(T(3.KE,I)-T(3.K.I))* 

A         T(  1 ,K  ,  I  ) 
I   NCY  =  ( Y-YMIN  )/0L Y 


lb     i\TEkSECTION     bITHiN     dOUNDS     Of     THE     SURFACE 

I  F(Y  .bT .T(2,^ .1)  .AND.Y  .bl.T C2,KE.  I) )  GO  TO 
I  F(Y  .LT.T(2. K.I)  .AND .Y  .LT.T(2.KE.l) )  bO  TO 
IF(ZcAR.LT.T(3.K.I).ASD.ZdAk.LT.T(3.K£,I)) 
IFfZL"AR.bT.T(3,^.n.AND.ZBArA.GT.T<3.KE.i'>) 
1  F(Y  .LE  .  AKP«A  X  (  I  t  )  )  bO  TO  1  V  :; 
A   K^'A  X  (  I  Z  )  =  Y 


1v5     i   F(Y  .bE  . AKMI  N(IZ  )) 

A  KMNMZ  )  =  V 
1V6     C  CNT  I\UE 

2  C    CONTINUE 

^5  ,C  CM  iNUc 

4"i     C   C.\T  lUUE 


bO     Tu     1 Vc 


d5 
25 
bO 

bU 


TO 
TO 


25 

25 


Lu      rANuof^     NCIbL      TG     AK>"--(li) 


^  L  0  1      h  1 1  r.  T  S      Y  L  ( I  ) 


I   FM  HLOT  .EQ.  .  ■»     uu     TO     c5 
Z   L(  I  )=ZL  AR 

Z   L  (  .  )  =  Z  ^  A  k 

h£N(AKrAX(in.LT.F'LT1.uw.ANf^Ax(IZ).bT.PLTL) 

LNLti>SrA<h  iNnZ">.LT.HLT1  .OK  .AM«»IN  HZ  )  .bT.PLTa> 
YL  (  1  )  =  /n^^  IN(  12  ) 
K  ^  I  N  =  K  f-  I  N  ♦  1 

CALL     CUKVtCtL.ZL.      1.-i) 
...FIN 
. .  f  i  i\  ' 

Lbc 

*HtN{AKvMlN(IZ).LT,PLT1.0n.AKMlN(IZ).bT.PLT2) 
VL  n  >=  A^fiAX  (  IZ  ) 

k;nax  =  k^ax*  1 

CALL     C  UhVE  (YL.ZL  .     1.-1) 
...FIN 
LLSE 

YL  (1  )=  AK«AX<  IZ ) 

Yl  (2)  =  AKMIN(  IZ ) 

Kl"  AX=<  f<A)t  ♦  1 
•        K^'IN=l(f"lN*1 
.        CALL     CLKVE(YL.ZL.     2  ,-  \) 
...FIN 


.   . .  F  I  r^ 

a^ 

C   ONT iNUt 

o5 

C  Cwl  INUC 

100 

C  UM  INUt 

L  LTR=LC  TR*1 

J  MN=0 

J   (<  A  X  =  0 

STAkT      iKDtX     FOR     Y-KASTcR 

L  l*A*  =0 

L  M  N  =  u 

w  f^ITE  (6  .  1222  )       LCTf< 

DC    i)CU     I     =1,lNCYh' 

Y  LArt  =  1     -I'LY     •♦  Yi-oIN 

UN  [J      Z     iNTEKSECTIoNb. 

d  Nf-A  x(I      )  =  -1..t1».^ 

e  Kf  I  N(  I      )  =  1U  .c  1C 

DO^Cu     J  =  1,NCi) 

I    F(X  NORf"  U  '>  .  L  t  .U  .  )     uO     TU     A  U  o 

DO    i(.U     K  =  1.A 

K  E  =  i^  +1 

I   f  (H  .ta  .  4)     <  t  =  1 
IF(Mli(T<c,K.J)-T(^,KE,J)).Lt.l.r-0)     00     TO     i^O 
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AKl<"--(i)  =  huKI20MAL     LIMtTS     fOh     SEub^^NTIAL     l      SCAN 
DKr>!--(J)=VtRTICAL     LIMITS     FoR     SuCiUbhTlAL     Y     SCAN 


CEMtK     Ii"AtE        yllh     RtSPECT  70        FlKST     iMAbE     ANu 

fUST    I^'     SlAN     Ih.MGE     bOTTOM  TO     TOP      (I)     ANU     LEFT 
DliC  KtT  fc"     I  .nTE  RV  ALS  . 

CONVLkT     iVA>.     ANb     MN.     EOVat  COO^UINATkS     TO     bKIu 

b   t  Y  I"  I  i\  =  Y  M  N  -  y  !»•  I  N  1 
0   tTi'iAXryrAX-YhAX"! 

DtLY  =  (    oEvf-I^•♦UEy^'AXW2. 

Y  UN  =Yr*INl-»D  ELY 

D  LYl=AcS(V!^AXl-YMiNl)/RNCYM 

0  cZi^lNsZf'IN-ZrtlNl 

0  tZMAXsZ^AX-ZK'AX  1  ' 

OELZ=(     DEZMI  N-»DEZf*AX  ) /2  . 

2  HIN  =  ZMIM*0  ELZ 

0  LZAl  =  AcS(ZMAx1-ZMlNl)/KNCYM 
I NCTh=0 

,  UO     105  0     I  =  1,lNCYt*! 

'  F   X=(  mk:«*A  X  (  I)  -YMIN)  /DL  Y  1 

1  X  =  I  F  I  X  (  F  X  ^ 

\  G  X=(  A^h1IN(  I)  -YMIN) /OL  Y  1 

!  I  b=I F IX  (6X  ) 

2  t*    CHARACTER  FLAG 
I  KF=C 

COKThAST  PIXEL  FLAb 

I  tF=: 

D  L    1 CoO     J  =  1,  INCYH 
F  Y=(bKf*AX  (  J  )  -2MN)  /DLZ  Al 
J   F=I F  IX  (  FY  ) 

G  Y=(bK«iN(J)-Zf^IN)/OLZAl 
J  ti=I  F  IX  (GY) 
I  R  =(  I  -1  )*  INC  Y?"*  J 

I  FdR.Nt  .0)     GO     TO     103  j 
I  >^  =  2A 
I  RF  =  1 
1  ciO    C  ONTINUe 

IF(IX.Nc.J)     GO     TO     1033 

CmH  D (IK  )  =  0N  b 
I  £F=1 

GOTO     1u44 
1  Li3     C  OnTINUE 

CAHD(IR)=ZEKO 
IF(IG.NE.J)     to     TC     10  30 

CAkD  f IK  >  =  0N t 
I  kF=1 

I  I,CTK  =  INCTR*  1 
G  0     r 0     1 LA4 
1  (.06     C  CNT  INUb 

CAkD  (  I  K  )  =  ZE  f«U 
IF(jF.Nt.I)     Cu     TO     IjAj 

CAkC (  I  K  )  =  0N  E 
I   cF='! 


USE     GHlu 
TO     H  IGH  T 

iNTt  Gt R  b  . 


INTERVALS 
(J)     AT 


OF 
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L    ^-(T(».».E.J  )-T(:.K.J))*(Ybnf<-Ti:.K,J))/(T(^.KF.J)-T(^,K>.J)) 
T(  i  .  Ik  ,  J  ) 


.i     IfjTthScCTION     wIlhlN     30UNU:i     0^     Tht     SUhFACt. 


c55 
4L0 


M  YL  A 
Ffi  f  . 

►..•   r,   X   ( 

j«  «  A  =  L 

t.-,T  IN 
F'4  F  . 
Kf  i  N  ( 
(•'  I  f,  =  L 
t.sT  IN 
UNT  IN 


LT  .T 
u  I  .T 
W  .  GT 
K  .LT 
L  t  .b 
I  )  = 
^.  /,  X  + 
Ut 

I    )  = 

Uc 
Ut 
UE 


(3. N, J). And. ZF, LT, T(  }.  ,tkt. J))     uC' 
'3  ,K  ,  J  >  .  AND  .ZF  .GT  ,T(  3  ,»  £,  J  >)     t.  0 

.T(i,K.J),AND.YBAr<.GT.T(2,KE.J)) 

.T(^.K,.J).AND.YBAh.LT,T(2.Kt,J)) 

K  >t  A  X  (  I      )  )     GO     TO     I  j^ 

2F 

1 


TO     30: 
TO    3  30 
b  0    T  D    i  G  u 

u 0   TO    3 n u 


K  «  I  fH  I 

ZF 

1 


)  )     GO     TO    c -  5 


ADO     KANuCf"^    NOISE     TO    3KM--(It) 


l-LOT      F-OInTS 


Fa 

LCt 
L(1 

HfcN 

U 


..F 

L:.t 


PLOT  . 

)  =  Yd  A 

)  =  Yu  A 

{  b  K  ?i  A 

lUE  SS 

ZL  ( 

J.'"  I 

C  AL 

.  .F  ir( 

IN 


EG . C)     00     TO     ^«5 

K 

X(I      )  .L T.PLT3.0R  .bN^AX  (I     )  .  G T  .  P L T A ) 
(bX^•lN(I     ).LT.PLT3.0R.eKMlN<I      ).GT.PlT(.> 
1  )  =  r  K.  M  N  (I     ) 

N  =  J  MN+  1 

l   coHvfcCYL.ZL.    *,-n 


or.  T 

Cr.  I 


rc  N 
ZL 
Jf^. 
CA 

.  .F  I 

LSE 
Zl 
Zl 
J  I" 
Ji^ 
Cm 

.  .F  I 

1 14 

If.Uc 

I.mUc 

I  i»OL 

^LC  r 


(eKMNd      ).Lr.PLT3.uR.cKMIS(I 

n  '>  =  bK^Ay  ( I    ) 
Ax  =  j  ^  mx-»  1 

LL     CUKVE(YL.ZL.      l.-l) 

N 

(l)  =  tK.t"AX(I     ) 

(  2  )  =  LKI-iIN(  I      ) 

AX=JfAX'»'» 

II>(  =  J  f'  I  N  ■»  1 

LL     CUKVECYL.ZL,      2.-1) 

N 


)  .bT.PLT^  ) 


♦  1 
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1  I  *.  3    C 


1.50 


I     -10 


iOC 


C    Tu     1 

LAr  D  f  I 
K  J  G  .  N 
CAnO  (  I 
tf  =  1 
0  TO  1 
UNT  INU 
CARD  (  I 
CNT  iNU 
CNT INU 
CNT  I.-,U 
CTK=LC 

h  1 14  T  2 
mLl  t  N 
CTk=LC 

gntinu 

CTh=LC 
ETUK.^ 

'ELEWA 
NO 


oA4 

t 

K  ^  =  Z  6  h  u 

t  .1  )     bO 

H  )  =  0N  £ 


TO     lUi 


)  =  ZE  f<U 


R*  1 

T  .EQ  ,c)     bO     TO     110 
C,  AZ  ,tL  t  PULL 
PL(U)  i 

R*1 


1H1. 

HON     ANGLE  =  '.  F3.^./  .  JX  , 


5X.'AZII«UTH    ANSLc= '.F8.2,/ .5X  . 
ROLL     ANbLE  =  ', F8,2  ,/  /  J 
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\JUO\Jtil     CONVEhTS     Ai^     iNf-UT      IhXN)     AhWAY      I'jTO     A     (^NX2n)     A  h  fv  a  y     bY     mDDING 
bACKGRi-UND     OATA. 

AkIN  =  INHljT     AHRhV,     Mk07=CUTPbT     Af<rt^V,     lND(i=i'»PUr     0  I  K  t  N  T  i  0  N 

UIMlNSION     ARIN(INUi>').     APGT(NOD) 
NO0t  =  <:«lNDM 
NOD  =i,*i  NUp^ 
FO^*  =  fLOAT(I^u^') 
ft  D  M  =  i  I.  R  T  (  F  0  ^  ) 
1DK  =  I  FI  X ( HOM  ) 
lOH  c  =  c*  I0i\ 
Z'J        FOH^  AT  (  ZX  .  '  I  uR     ='.,  15./) 

j  =  N(jw  iNDbx  FOR  Output  array   i  =  culunn  ini)Ex 

DO  J,uU  J  =  1.IUR£: 
I  F ( J  ,GT  .IDR  )  GO  TO  270 
DC  ZuQ  1=1. IDR 
I  J=  Il;R2*(  J-1  )♦! 
I  JK  =  IDR*(  J-1  )  •»! 
AP0T(IJ)=ARIN(IJK) 
2 ..  1   C  0  N  TI  N  U  £ 


FILL  REMAINDER  uF  ROw  WITH  ZEROS 

c2<,     DO  Zbu  1  =  1,  IDR 

I  J=  ll;Rii*(  J-1  )-»iOR-»  I 

A  R  0  T  (  IJ  )  =  0  . 
25  u    CONTINUE 

GO     TU    3  GO 
27U     CCNTINUfc 

FILL     Rfci^AINDtR     OF     RO-3     yITH     ZEROS 

DC     lab    1=1. IDR.; 

I  J=  IUR2*( J-1 ) ♦ i 

A  R  0  T  (  I  J  )  =  0  . 
<:eo  CON  TINUE 
iOo  CCNTlNUb 

RETURN 

END 
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THIS  SUBROUTINE  CALCULATES  THE  CORRELATION  BETWEEN 

THE  A  VECTORCOFT  C OE f f  1 C  I  ENTS  )  AND  THE  B  VECTuR^DFT  COEFFICIENTS) 


2037 
63 


103 

10^4 

20o 


71 
72 
73 


SUBROUTINE   CALCOR(KK2.JJ2.NV.AZ.tL.ROLL,J  .N12.N22 

DIMENSION  yH (16384)  ,WHR(16384  ) 

DIMENSION  5(128), INV(12a),Mf3) 

DIMENSION  CMAX(20) 

DIMENSION  MAXX (20)  .MAX  V  (20) 

DIMENSION  W(  12ti,12£>) 

C0MPLEX*8  A f 16id4)  ,e( 16384) 

COMMON  WHR.UH.A  ,CMAX 

M(1  )  =  KK2 

M(2)=JJ2  t 

MC3)=0 

F0RMAT(/.5X, 'J  =',I2,3X,'CMAX  =',F17.6./) 

F0RMAT(5X,'n12  =',IA,'N22  =',I4,'NV  =',I4) 


) 


N V  =  INPUT  D  IMENSIO  N 


DOUBLE  SIZE  OF  DATA  ARRAY 

yHR  =  INPUT  VECTOR.  WH   =  OUTPUT  VECTOR 

NV4  =4*NV 

CALL  VDOUBL  (kiHN  ,WH  ,NV  ,NV4) 


PLOT  DATA  ARRAY 

TRANSFER  THE  WH  <K)  DATA  TO  THE  COMPLEX  B  MATRIX 

NDATA=0 

DO  104  l2==1fN22 
DO  103  I1=KN12 
IWH=     (I2-1)*N12*I1 

B (I WH)=UH     (I WH) 

IF(     WHMWH)      .GT.0.0)     NOATA  =  ND  ATA*  1 

CON  TINUE 

lWHM=lWH-9 

WRITE(6,72)   (WH  (I)  ,  I  =  IWHM. IWH) 

CONTINUE 

WRITE(6.208)  NOATA 

FORt"AT('   NDATA  ='.17) 

I F( NDATA.EQ .0)  STOP 

COMPUTE  THE  DISCRETE  FOURIER  TRANSFORM  OF   THE  wH  (K)  DATA 

FORMAT( 16(1X , F7.2) ) 

FORI-AT  (10(1X.E12.5)) 

FORMAT(7(2X.  15)  ) 

CALL  HARM(B.M,INV,S  .1  .IPERR ) 

COMPLEX  MULTIPLY 
DO  22U  11=1, NV4 


on 


220 


e(ii)=A(n)*coNj6(b(ii)) 

IFCn  .LT.N22)  WRITE(6,101) 

CON  TINUE 


3Q 


40 


41 


II  .A(ii),e( II ) 


INVERSE  Fourier  tkansforw 


CALL  HARM(9,M,INV.S ,-1 .IPERR) 


CALCULATE  THE  MAGNITUDE  OF  THE  FOURIER  COHPONENTS 
AND  STORE  ONLY  THE  POS  FREQ  VALUES  OF  d(I1,l2) 


DO  30  I  2=1, N22 

DO  3U  11=1. N12 

IWH=  (I  2-1>*Nl^Ml 

w<I  1f I2)  =  CABS  O  (IWH  )) 

IF( W(I1 .12) .LT.CMAX  (J  ))  60 

CMAX(J)=U(I1 .12)   . 

MAX  X(J )=I  1 

MAX  Y<J )=I2 

CONTINUE 


TO  30 


GO  TO  41 


42 


PRINT  THE  OFT  COEFFICIENTS 

IF<MAXX(J ).6T.3)  60  TO  40 
MM3=1 

MM2  =  2 

MM1=3 

MAX  X (J) =4 

MP1=5 

MP2=6 

MP3  =  7 

MP4=b 

GO  TO  42 

CONTINUE 

Nl2M=Nl2-4 

IF(MAXX (J ).6T.N12M) 

MM3=MAXX( J)-3 

MM2=MAXX( J)-2 

MM1 =MAXX< J)-l 

MP1  =MAXX(  J)-»1 

MP2=MAXX(  J)-»2 

MP3=MAXX( J)*3 

MP4=MAXXO)*4 

GO  TO  42 

CONTINUE 

MM3=Nl2-7 

MM2=Nl2-6 

MMI  =M2-5 

MAX  X( J)=Nl2-4 

MPl=Nl2-3 

MP2=Nl2-2 

MP3=N12-1 

MP4=M2 

CONTINUE 
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10U 
10  1 


61 


kRI  Tt(6.1Ca)  Mn3.Mtt2•^'M1.^»AXX  cJ),f«Pl,MP2fMH3f^PA 

FCRrATC       K',5X.'A('.I3.'.K)',7(6X,'A(',I3.'.K)')) 

f  ORf^AT<2X,I4,8(lX,t12.5)  ) 

DO  61  12=1  ,Ni2 

ttRITE(6.1U1)l2,(W(I1.I2).I1=MM3,MP4) 

CON  TINUE 

RETURN 

END 
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KAL  IS  A  RADAR  RANGE  ESTIMATOR  USING  A  KALIAN  FILTER. 
THERE  ARE  2     RANGE  MEASUREMENTS  AND  A  RADIAL  ACCELFRATION 
MEASUREMENT. 


ESTIMATION 


MODEL 

x(n*1 ) 

Y(K*1) 

KALMAN 
X'(K*1 
X-HAT( 

P(K*1) 
P(C)  = 
T  =  TR 
SUBROU 


DIM 

*TX( 

DI»« 

*RNO 

DIM 

DIM 

DIM 

DIM 

DIM 

COM 

•  DEL 

COM 

♦ROL 

*RAN 

•ANA 

OAT 

OAT 

DAT 

DAT 

DAT 

1/(S** 

CAL 

Nl  = 

N0  = 

TGS 

TGS 

DT  = 

NP  = 

NPM 

NPH 

RES 

R90 

R27 

6YD 

NRU 

10  CON 


=  H(K)*X 

=  S<K*1) 

ESTIMATE 

)  =  H(K)* 

»C*1)  =  X' 

=  P'(K*1 

)  =  H<K)* 

=  P'(K*1 

G(0)*Q(0 

ANSPOSE 

TINES  USE 

ENSION  X( 
4).Y(3),R 

ENSION  TP 
152(100) , 
ENSION  AP 
ENSION  WH 
ENSION  PA 
ENSION  10 
ENSION  IP 
MON  XT,YT 
T,NPRF,TR 
MON/lNTCy 
406  ,T6S1 
GTM,XR,YR 
,ANE,PAUV 
A  IN,N,IT 
A  IS, Ml/2 

A  IL.L/2* 
A  ISEED/A 

A  PI/3. 14 
2)  "ODEL 
L  FR80ID( 
5 
6 

1  =  10. 

2  =  10. 
0.05 
31 

=  NP-1 
=NP/10 
=  0.5 
=PI/2. 
0=3.*PI/2 
=10.7252 
N  =  0 
TINUE 


(K)*6(K)*U(K) 
♦X(K*1 )*V(K*1 ) 

X-HAT(K) 

(K)-K(K*1)»(S(K-H)*X'<K)-Y(H*1)) 

)*ST(R*1)*{S(K+1)«P'<K*1)*ST(H*1)+R(K*1) )• 

P(K)*HT(K)+G(K)*Q(K)*6T(K) 

)-K (K*1 )*S(K*1)*P'<K*1) 

)*GT(0) 


-1 


0:    INTCON 

MANVR 

F7KALI 
4),H(4,4),6(4),S(3»4)  ,P(4,4),T1(4,4),T2(4,4),T3(3), 
(3,3),F6(4,4) 

(100) ,EP1(100),XP1<100) ,RN0IS(100),VP1 (100) ,VP2(100), 
XTP(100>,2TP(100) 
1(100),AP2(10O),AP3(1OO) 
EN(2)  ,TME0D(2) 

UV(3),ANUV(3),ANL(3),ANLR(3),TPA(3) 
PT(5),STAT(5) 


/3*4/ 
*3/ 
1/ 
2573/ 

1592654/ 

STATE  VECTOR 

'  6220  R.H.  WILLIAMS 


6262  ') 


3-D  TARGET  TRAJECTORY  DATA 

ORIGIN  AT  INITIAL  TARGET  POSITION 
RT  =  TRUE  RANGE  RELATIVE  TO  RADAR 
RDOTT  '     TRUE  VELOCITY  RELATIVE  TO  RADAR 
»?DOTT  =  TRUE  ACCELERATION  RELATIVE  TO  RADAR 
VELOC  =  TARGET  VELOCITY  RELATIVE  TO 

TRAD-  =  TURN  RADIUS  OF  TGS-  6. 

TRAT-  =  Turn  RATE     OF  TGS-  6.  TURN    (PAD/SEC) 

-R  =  RADAR  LOCATION  AT  T(0) 

-T  =  TRUE  X,V,2  TARGET  POSITION 

V-T  -  TRUE  X,Y,Z  VELOCITY 

A-T  =  TRUE  X,Y,Z  ACCELERATION 

A-T6  =  TRUE  X,Y»Z  ACCELERATION  (G'S) 
=  VERTICAL  (A7)  ACCELERATION  RELATIVE 


ORIGIN 

TURN 

TURN 


ANA 

ANE  =  HORIZONTAL  (£L)  ACCELERATION 

PAUV(I)  =  PITCH  AXIS  UNIT  VECTOR 

PAUV( 1 )  =  X  i 

PAUV(2)  =  Z 

PAUV (3)  =  Y 


TO  TARGET  OPIENTATION 
RELATIVE  TO  TARGET  ORIENTATION 


ATC  = 
TCK  = 
ATC 
TC 
TCK 
Q  =  CO 
=  EC 
REA 
DO 
1F( 
DO 
00 
H  =  ST 
DAT 
A6»«TK 
ART 

F( 
F( 
F  ( 
F  ( 
F< 
F( 
H(1 
H(1 
H(1 
H(1 
H(2 
H(2 
H(2 
H(3 
ANLRK 
ANL 
H(4 
G  =  ST 


TI«E 
TRAN 

=  2 
'     AT 

=  ( 
VARI 
U*UT 
0(5, 
115 
lOS. 
118 
116 
ATE 
A  H/ 
=  6A 

NOR 
IH.E 
IH.E 
IH.E 
IH.E 
IH.E 
IH.E 
•  1)  = 
»?)  = 
,3)  = 
,A)  = 
,2)  = 
,3)  = 
.4)  = 
.3)  = 
=  RA 
RK=1 

.^)  = 

ATE 


CON 
SITI 

• 

C/DT 

1.-E 

ANCE 

3  = 

150) 

I0S  = 

EQ.2 

IR  =  1 

IHs:1 

TRAN 

16*0 

««^A 

«AL 

Q.I) 

0.2) 

a. 3) 

Q.I) 
Q.2) 

U.3) 

1. 

DT 

A6MT 

ARTK 

1  . 

DT»A 

ARTK 

1  . 

DIAL 


STANT  OF  ACCELERATION 

ON  MATRIX  ACCELERATION  TI»E  CO^STANT 


XP( 

MA 

(RM 

Q» 

1.2 

)  D 

.3 

.2 

SIT 

./ 

ACC 

RAD 

AG 

AG 

AG 

AR 

AR 

AR 


-TO) 

TRIX  OF  U  (LXL)  -  INPUT 

S  RAN6E/SEC*SEC*StC )**? 

RCV.ILAST 

T  =  0.1 


ION  MATRIX  (NXN) 

ELERATION  TRANSITION  MATRIX  CONSTANT 

lAL  LOAD  ACCEL.  TRANS.  MATRIX  CONSTANT 

MTr  =  1  .0 

MTKsO.O 

MTt(=0.5 

TK=0.0 

TK=1.0 

Tk=0.5 


K*0T*DT/2. 
*DT*DT/2 . 

QWTK 
•DT 

NORMAL  LOAD  ACCELERATION  TRANSITION  CONSTANT 


ANLPK 

NOISE  INPUT  MATRIX  <NXL) 


34 


RDOTT  = 


D 

I 
I 
I 

P 

R 

R 

D 

MANV 

VYT, 

LOC 

*     = 

Y  = 
Z  = 

C 

V  = 

V 
A 
A 
A 

AT  = 
A 
A 
A 
A 
A 
A 

AH  = 
A 
F 
E 

VH  = 
V 
A 
A 
E 
E 
A 

1C7S 
I 
R 


0  100  11=1 

P=I1*1 

PS  =  0 

F(II.EQ.I) 
CTNS=0.7 
NOIS(II)=S 
N0IS2(II)= 
ELR=WELOC* 
R  CALCULAT 
VZT)  TAR6E 
ATION  RT  0 
DOWN  RANGE 

UP   RAN6E 

OFF  RANGE 
ALL  MANVR 
TOTAL  TRUE 
=SQRT (VXT* 
XT=(VXT-VX 
YTs(VYT-VY 
2T=(VZT-VZ 

TOTAL  TRU 
T=SQRT(AXT 
XT6=AXT/10 
YT6=AYT/10 
2T6-AZT/10 
NE*((AYT*V 
NA=(AZT*VZ 

HORIZONTA 
HsSQRT (AXT 
LE=ANE*COS 
LN=SQRT (AN 

HORIZONTA 
H=SQRT(VXT 
2=ATAN?(VZ 
Z0E6=AZ*18 
L=ATAN2(VY 
L0EG=EL*1S 
NAEL=ANA/E 

=  COOROIN 
F  <  ICTS.EQ. 
0L5  =ATAN? 


,NPW 


IPS  =  1 


QRT(R(1,1))*66N0F(1SEED)*PCTNS 

SQRT(R(2,2))*6  6N0F(ISEE0)*PCTNS 

0ELT*FL0AT(NPRF) 

ES  TRUE  POSITIONS  (XT,YT,Z7)  AND 


TRUE  VELOCITIES  (VXT, 


STARTS  3  ORIGIN  S  T'O,  AND  FLIES  TOWARD  RADAR  FIXED 

X  AXIS 

REFERENCED  TO  TARGET  3  T  =0 

REFERENCED  TO  TARGET  S  T  =0 

REFERENCED  TO  TARGET  3  T  =0 


COORD  SYSTEM* 


T(0) 


TO  COORD  SYSTEM 


T(0) 


VELOCITY  RELATIVE  TO 
VXT+VYT*VYT*VZT*VZT) 
TM)/(DELT*FLOAT(NPRF)) 
TM)/(DELT*FLOAT(NPRF)) 
TM)/(DELT*FLOAT(NPRF)) 
E  ACCELERATION  RELATIVE 
•AXT*AYT«AYT*AZT*A2T) 
.72S2 
.7252 
.7252 
XT-VYT*AXT)*VXT-(AZT*VYT-VZT*AYT)*VZT)/(V«VH) 

T-VZT*AXT)/VH 

L  ACCELERATION  -  RELATIVE  TO  COORD.  SYS.,  T(0) 

*  A  X  T  ♦  A  Z  T  ♦  A  Z  T  ) 

(EL) 

A*ANA*£LE*ELE) 

L  VELOCITY 

*VXT*V2T*VZT) 

T.VXT) 

O./PI 

T,VH) 

O./PI 

LE 

ATED  TURN  SWITCH 

1)  ROL4=ATAN2(ANA,10.7252) 

(ANA, 10. 7252) 


-  RELATIVE  TO  COORD.  SYS.,  T(0) 
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G  (  1  )  =  DT*DT*OT/A, 

G(2)=DT*DT/2, 

G(3)=0T 

6(4)=     6(3) 
R  =  COVAHIANCF  MATRIX  OF  V  ( «1 1  X  f*  1  )  -  i«  E  A  S  UR  E  ^»E  N  T 

DATA  R/9*0./ 
R(1,1)=RCV 
R(2,2)=      R(1,1) 

R(?.3)=100*R(1,1) 
S  =  OBSERVER  MATRIX  ( « 1 X N )    Y  =  SX 

DATA  S/1?*0./ 

S(1  ,1)  =  !. 

S(2,1)=1. 

S(3,4)=1. 
INTCON  INITILIZFS  TRAJECTORY  VARIABLES 


THRESHOLDS 


CALL  INTCON 

TX  =  ITERATED  EXTENDED 

TX(1)  =  1  . 

TX(2)=TX(1) 

TX(3)*TX(2) 

TX(4)sTX(3) 

K=0 

T  =  0. 

LL*0 

INCP=0 

RN0IS(1)*SQRT(R(1,1))*66N0F(ISE£D) 

RN0IS2(  1)=SQRT(R(2,2))*66N0F(ISEED) 
y  =  INPUT  OBSERVATION  VECTOR  (MlXD 
MULTIPLE  MEASUREMENTS 

Yd  )  =  RT  +  RN0IS(1) 

Y(2)=RT*RNOIS2(1) 

Y(3)*0. 
PRESET  ESTIMATOR  INITIAL  STATE 
X  (1)  ^     RANGE (K) 
X(2)  =  RANGE  RATE(K) 
X(3)  =  RANGE  ACCELERATION 
X(4)  =  ESTIMATED  RADIAL  NORMAL  LOAD  ACCELERATION 

X100FF=  1. 

X20k=1.1 

X  (1  )  =  Y(1 )*X100FF 

X(2>=X20K*RD0TT 

X(3)=0. 

X(4>=0. 

ePl (  1)=  RT-X(1) 

XP1(1)=X(1) 

VP1(1)*X(2) 

VP2(1 )=RDOTT 

AP1 (1)=x(3) 

AP2(1)=0. 

AP3(1)=0. 

EPMAX=0, 

RDRX=RT 

NR(JN~NRUN'^1 

WRITE(6.300)  Q,RCV,RT,R0OTT,A6MTK,ART(t.S(3.3),NRUN 


1A4  FOR"<AT('   COOD.  TURN  SU  ='I3,'   BOH  =',F9.3,'   COORD.  TURN  ROLL 
*',F9.3) 
R0LDEG  =  R0L4*1iJ0./Pl 

ANADOT=((-VXT«AXT-VZT*AZT)*ANA)/(VH«VH) 

AZOOT=ANA/VH 

ELDOT=ANE/V 

R0LD0T  =  10, 7252*  AN  ADOT/dO. 7252*10. 7252-»  AN  A  *ANA) 

V-TM  =  DELAYED  TRUE  X,Y,2  VELOCITY 
VXT«sVXT 
VYTM=VVT 
VZTW=VZT 

-3  =  X.Y.Z  DIFFERENCE  BETwEEV  RADAR  AND  TRUE  TARGET  POSITION 
X3=XR-XT 
Y3=YR-YT 
Z3=ZR-ZT 

RAN6T  =  TRUE  RANGE   RANGTM  =  DELAYED  TRUE  RANGE 
RAN6TM=RAN6T 

RAN6T=SQRT(X3*X3+.Y3*Y3-*Z3*Z3) 

VNOIS=RNOIS (II) 

RAN6N=RAN6T+VN0IS 

RAN6N2=RANGT+RN0IS2(II) 

RE<4DR=AMOD(RAN6N,RES) 

REi^DR2  =  A^0D(RAN6N2»RES) 

RES2=RES/2. 

RANGE  -  MEASURED  RANGE 

RANGE=RANGN-REnDR 

RAN6E2  =  RAN6N2-R£?^DR2 

RT=RAN6T 
AWNL  =  MAGNITUDE  OF  NORMAL 

AL=13. 

BE*-3, 

6AM=0,5 

EPS=6GN0F (ISFED) 

AMNL*AL*BE*EXP(6A«*EPS) 
AZ  '  AZDEG 
EL  =  -EL0E6 
ROLL  *  -R0LDE6 
PAUVd)  *  X 
PAUV(2)  =  Z 
PAUV(3)  =  Y 

SAZ=S1N(  AZ) 

CAZ=C0S(   AZ  ) 

SEL=SIN(-EL) 

CEL=COS(-  EL  ) 

SR0=SIN(-R0l4) 

CR0=C0S(-R0L4) 

TPA(1)=PAUV(1)*CEL*CAZ-PAUV(2)*CEL*SAZ-»PAUV(3)*SEL 

TPA(2)=PAUV(1)*(CRO*SAZ*SR0*SEL*CAZ)*PAUV(?)*(CR0*CAZ-SR0*SEL*Sa; 

.-PAUV(3)*SR0*CEL 
TPA(3)=PAuV(1)*(SRO*SAZ-CRO*SEL*CAZ)*PAUV(2)*(CRO*SEL*SA2*SRO*CA: 

.♦PAUV(3)*CR0*CEL 
ANUV(I)  =  NORMAL  ACCELERATION  UNIT  VECTOR 

=  (TPA(I)  X  V-T)/V     (VECTOR  CROSS  PRODUCT) 
(X,Z  ,Y) 
ANUV(1)-((TPA(2)*VYT)  -  (TPA<3)*VZT))/V 


WITH  NOISE  AND  QUANTIZED  TO  RESOLUTION 


LOAD  ACCELERATION  (  G'S  ) 
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CN  THE  RADIAL 


ANUV(2)  =  ((TPA(3)*VXT)  -  (  TP A  (  1 ) • V YT ) ) / V 

ANUV(3)=((TPA(1)*V2T)  -  (TPA(2)*VXT))/V 
C    ANL(I)  =  NORMAL  LOAD  ACCELERATION 

ANL(I)  =  ANUV(1)*AWNL 

ANL(2)  =  ANUV(2)*Af*NL 

AML(3)  =  ANUV(?)*AWNL 
C   ANLR(I)  =  RADIAL  NORMAL  LOAD  A C C £ L E R AT  I  ON (  6'S  ) 
C  =  PROJECTION  (DOT  SCALAR  PRODUCT)  OF  ANLCl) 

C  VECTOR  FROM  TARGET  TO  RADAR/RADIAL  RANGf 

ANLR(1)  =  ((ANL(1)*X3)>(ANL(2)*Z3)*(ANL(3)*Y3))/RA**6T 

Y (1 )=RAN6E 

Y (?)=RANGE? 

Y (3)=ANLR  (1  )*6YD 
C    VR2D  =  1ST  DIFFERANCE.  RANGE  VELOCITY   VP200  =  DELAYED  VR2D 

VR20D=VR2D 

VR2D=(RAN6T-RAN6TM)/(DELT*FL0AT(NPRF)) 
C   AR  =  RANGE   ACCELERATION,  2ND  DIFFERANCE 

AR«=AR 

AR=(VR20-VR2DD)/(0ELT*FLOAT(NPRF)) 

WRITE(6,133)  RAN6T,RAN6T«,VR2DD,VR20,AR,Y(3),X(3), 
*X3,Y3,Z3 
C   TEST  RADIAL  2ND  DIFFERENCE:  AR,  ft  LIMIT  TO  »«AX,  EXPECTED 
C   RATE  OF  CHANGE:  SQRT(Q) 

ADI F=ABS(AR-ARW) 

OLI«=SQRT(Q)*DT 

IF(ADIF.6T.QLI«)   AR  =  S16N(QLIM,AR)  +  AR^' 
C    TEST  FOR  RE-INITILIZ ATION  OF  KAL«AN  FILTER 

S16THS2.0 

S161  =  SI6TH*SQRT  (R  (1  ,1  )) 

SI62=SI6TH*SORT(R(2,2)) 

SI63  =  AMAX1 (SI61  ,SI62) 

PHAXsQ* 

EEST=((V(1)*Y(2))/2.)-X(1) 

IF(EEST     .LT,SI63)  60  TO  44 

INCP=INCP*1 

IPS«1 

00  42  JP=1,N 

PMAXSAMAX1 (P(JP,JP)  ,P«AX) 

CONTINUE 

PPCT=1.0 

DO  43  JP=1,N 

P  (  J  P  ,  J  P  )  =  F»  M  A  X  *  P  P  C  T 

CONTINUE 

CONTINUE 

WRITE (6,133)  AZ,EL,ROL4,TPA(1),TPA(2).TPA(3),ANUV(1),ANUV(2). 

.ANUV(3),A«NL,ANLR(1)fAR 

CALL  FTKALI (K  ,X,H,6,Y,S,0,R    ,P , I N  ,  I S  ,  I L ,N ,«1 , L , T 1 , T 2 , I T , T3 , 
•TX,F6,IER,IPS) 

TsT*DT 

EP1(IP)=  RT-X(1) 

XP1  (IP)=X(1  ) 

VP1 (IP)=X(2) 

VP2( IP)=VR2D 

API (IP)=X(3) 

AP2(IP)=AR 


42 


43 
44 
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134 


A 
T 
X 
7 
V 
A 
P 
W 
F 
I 
I 
K 
I 
W 
I 

w 
u 

D 
F 
C 

W 

*v 
100  c 


194 
51 


P3  (I 
P  (  IP 
TP(I 
TP  (I 
PCT  = 
PCT  = 
CTW  = 
RITE 
ORWA 
F  (AB 
F  (AB 
=  K*1 
F(F6 
NS  =  F 
F(k»N 
N-SQ 
NHZ  = 
R=  (F 
ORMA 
ONTI 
RITE 
R?D 
ONTI 


P)=Y  ( 

)  =  T 

P)=)(T 

P)  =  ZT 

X  (2)* 

X(3)* 

T?(1  , 

(6,13 

TdX, 

S(EP1 

S(EP1 


3) 


0T*100./EP1(IP) 
DT*0T*50,/EP1 (IP) 
1  )*100./EP1 (IP) 


1  )*100./EP1 (IP) 

'4)  II,VPCT,APCT,PCTM,T2(1,1),PMAX 
13,'  PCT  VFLOV.  ACC,  («EAS  =  '  ,  1  0  (  F  1  ?  .  3  )  ) 
(IP)).6T.ABS(EPMAX))  EP«AX=EP1(IP) 

I  (1P)).6T»40.)  LL  =  LL*1 


(1  .1  )  .GE.1  .)  60  TO  51 

G(?,1)/( (1.-F6(1 ,1))*DT) 

S.LT.O.)  60  TO  51 

RT  (MNS  ) 

taN/?.*PI 

6(1,1)-DT*FQ(2,1))/(2*(1.-F6(1,1))*0T*WN) 

TC   NATURAL  FREQ.  (HZ)  =',F9.2,'   DA'^PING  RATIO  =',F8,3) 

NUE 

(6,130)  K,RT,y(1),X(1),')((2).EPl(IP),RN0IS<II),AXT6,AyT6,A2T< 

NUE 


I 
I 

U 

c 

w 
102  F 

w 

104  F 

I 

I 
C 

c 
c 
c 


0PT(1)=1 

0PT(5)=1 

PP=NP*1 

ALL  BEIU6R 

RITE(6,102 

ORMAT( 

;   VARIANC 

;   WAX IRUW 

R1TE(6,104 

ORl«AT(  /  ,' 

NO.  TIMES 
RPL0T=1 

F( IRPLOT.EQ.O)  60  TO  113 
ALL  INTAXS 
ALL  yAXAN6(0.) 
ALL  SIMPLX 
ALL  TITLE( 


(EP1  ,NPP,IOPT.STAT,IER) 
)  STATd  )  ,STAT(5)  ,EP«AX 

ARITHMETIC  MEAN  OF  RANGE  ERR  =',f7.3, 
E  (AVE.  RMS  ERROR)  OF  RANGE  ERROR  =',F9. 

RANGE  ERR  =',F9.3,/) 
)  LL,INCP 
0.  OF  Rl 

P  INCR.  =',I4) 


100, 'TIME 
ALL  XTICKS 
ALL  YTICKS 
0R6=AMIN1 ( 

MAX6=AMAX1 
ALL  GRAF (X 


ALL 
ALL 
ALL 
ALL 
ALL 
ALL 


CURVE ( 
CURVE  ( 
MESSA6 
RE  ALNO 
MESSAG 
RE ALNO 


IN 

(5) 

(5) 

ARM 
(  AR 

ORG 
ARM 

£p^ 

RNO 

cc 

( 

(' 

(R( 


KALMAN  FILTERS', 100, 
'ERROR  IN  YAROSJ', 
SEC. S', 100, 6. ,8.) 


IN(EP1,NP),ARMIN(RN0IS,NP)) 
MAX(EP1,NP),ARMAX(RN0IS,NP)) 
, 'SCALE', XMAX6, 

IN(TP    ,NP)  , 'SCALE', AR«*AX(TP   ,NP)) 
,   TP,NP,10) 
IS,TP,NP,10) 
OV  Q  =$',100,1  .,  b.9) 
Q,0, 'ABUT', 'ABUT') 
COV  R  *J',100,     'ABUT', 'ABUT') 
1  .1)  ,2, 'ABUT', 'ABUT') 


P.Q 
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114 


ALL 

ALL 

ALL 

ALL 

ALL 

ALL 

DUWM 

ALL 

ALL 

ALL 

ALL 

ONTI 

vplo 

F  (IR 

F(IV 

ALL 

ALL 

ALL 

ALL 


"ESSAG 

INTNO( 

DATE  (W 

TOFOAY 

MESSA6 

HESSAG 

T=LlNe 

LINESC 

LINES( 

LEGEND 

ENDPL( 

SUE 

T  =  0 

•  EQ.3 
PLOT.F 
INTAXS 
YAXAN6 
SIMPLX 
TITLE( 


('  # 
NRUN 
HEN) 
(TWF 
(WHE 
(T«E 
ST(I 
'ERR 
'NOI 
(IPA 
0) 


=i  .      100, 
.'ABUT', 'ABUT') 


00) 

N»10,1 .,9.3) 

00, 10, 'ABUT'  , 

PAK,150,bO) 

0RJ',IPA»C,1) 

SE J',IPAK,2) 

K    ,2,4.0,2.0) 


APUT','AbUT') 


ABUT') 


)  1VPL0T=1 
Q.O)  60  TO  114 


(0.) 


100,'TI»£ 
ALL  XTICKS 
ALL  YTICKS 
0R6  =  A«IN1  < 
MAX6=AWAX1 
ALL  GRAF(X 


ALL 
ALL 
ALL 
ALL 
ALL 
ALL 
ALL 
ALL 
ALL 
ALL 
ALL 
ALL 


CU 
CU 
ME 
RE 
ME 
RE 
ME 
IN 
DA 
TO 
ME 
ME 


OU«MY= 
ALL  LI 
ALL  LI 
ALL  LE 
ALL  EN 
ONTINU 
APL0T= 
F (lAPL 
ALL  IN 
ALL  YA 
ALL  SI 
ALL  TI 


RVE  ( 

RVE  ( 

SSAG 

ALNO 

SSAG 

ALNO 

SSAG 

TN0( 

TE  (W 

FOAY 

SSAG 

SSAG 

LINE 

NES( 

NES( 

6EN0 

0PL( 

E 

1 

OT.E 

TAXS 

XANG 

MPLY 

TLE( 


IN  S 

(5) 

(5) 

ARMI 
(ARM 
ORG, 
ARMI 
VP1, 

('CO 

( 

( 

(R(1 

('  # 

NRUN 

HEN) 

(T«F 

(yHE 

(TME 
ST(I 
'X(? 

'VR2 
(  IPA 

0) 


KALMAN  F  ILTER J' ,100, 
'   VELOC  IN  Y.ARDS/SECS', 
EC. $',100, 6. ,8.) 

N(VP1  ,NP)  ♦ARMIN(   VP2»NP)) 
AX(VP1 ,NP)  ,ARMAX(    VP?,NP)) 
'SCALE', XMAX6, 
N(TP    ,NP)  , 'SCALE', ARMAX (TP 

TP,NP,10) 
TP,NP,10) 


,NP)) 


»TP,NP,10) 

)V  0  =$',100,1.,  6.9) 
0,0, 'ABUT', 'ABUT') 
COV  R  =$',100,     'ABUT', 'ABUT') 
1,1), 2, 'ABUT', 'ABUT') 
t     =$',       100,     'ABUT', 'ABUT') 
J, 'ABUT', 'ABUT') 

00) 

N,10,1 .,9.3) 

00, 10, 'ABUT', 'ABUT') 

PAX, 150, bO) 

)$',  IPAK,1) 

0$',IPAK  ,2) 

K     ,2,5.0,1  .0) 


100,'TI«E 
ALL  XTICkS 
ALL  YTICKS 


Q.O)  60  TO  114a 
(0.) 

KALMAN  FILTERS' ,1C0, 
'  ACCEL  IN  YDS/SEC 'SECJ', 
IN  S£C.$',100,6.  ,?.) 
(5) 
(5) 


X0(J6  =  AMIN1(ARM1N(AP1,NP).ARMIN(    AP^.NP)) 
XWA)(6  =  A-AX1{ARWAMAP1,NP),AR»«AX(    *P?,NP)) 
CALL  6RAf(X0R6. 'SCALE'. XMAXG, 

AR«IN(TP    ,NP)  , 'SCALE', AR^'AX  (TP 
CURVECAPI,   TP,NP,10) 
CURV£(AP2,TP,MP,1C) 
CURVE(AP3,TP,NP,10) 
MFSSAGCCOV  Q  =S', 100,1.,  a.V) 
REALNOC      Q,0.'APUT','ABUT') 
MESSAGC   COV  R  =$',100,     'AeuT','ABUT 
REALN0(R(1,1),2,'AeUT'.'ABUT') 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


,NP)  ) 


'ABljT','AeUT') 


MESSAGC     #     =S'.  100, 

INTN0(NRUN, 'ABUT', 'ABUT') 
DATE (WHEN) 
T  0  F  D  A  V  (  T  ("  E  0  D  ) 
M£SSAG(WHEN,10,1 .,9.3) 
MESSA6{T^E0D, 10, 'ABUT', 'ABUT') 

IDUWI*Y  =  LINEST(IPAK,150,80) 

CALL     LINES('X(3)$',IPAK .1) 

CALL  LINESCAR   S',IPAK,?) 

CALL  LINES('V(3)$',IPAK,3) 

CALL  LEGENDdPAK  ,3,4.0,?.0) 

CALL  ENOPLCO) 
1U8  CONTINUE 

CALL  INTAXS 

CALL  YAXANGCO.) 

CALL  SIMPLX 

CALL  TITLE ('TARGET  GROUND  PLOT     l',100, 
•'  OFF  RANGE  (Z)  -  YDS. S', 100, 
*'  DOWN  RANGE  (X)  -  Y DS . $ ' , 1 00 , 6 . , 8, ) 

CALL  XTICKS(5) 

CALL  YT1CKS(5) 

XORG=        AR«IN(ZTP,NP) 

XWAX6=        ARWAX (ZTP,NP) 

YHAX6=AR«AX(XTP,NP) 

CALL  GRAF (X0R6, 'SCALE', XMAX6, 
*  ARmN(XTP  ,NP)  , 'SCALE' .YMAXG) 

NPLS=  5 

CALL  CURVE (ZTP,XTP,NP,NPLS) 
DRAW  ♦  AT  RADAR  LOCATION,  XT=RT  ,  ZT=0. 

CALL  MARKER (3) 

CALL  SCLPIC(2.) 

RDR Y=0. 

I F  (RDRX  .LE  .  YMAX6)  CALL 

CALL  RESET('SCLPIC') 
RESET( 'MARKER') 

MESSAGCCOV  Q  =$',100,1.,  fc.9) 
REALNO(      0,0,'ARUT','AfaUT') 
MESSAGC   COV  R  =S',100,     'ABUT 
REALN0(R(1,1),2,'ABUT','ABUT') 


CURVE  (RDRY  ,RDRX  ,1,1) 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


I 


ABUT') 


MESSAGC  *  =J',       100, 

INTN0(NRUN, 'ABUT', 'ABUT') 

DATE(WHEN) 

TOFDAY (T«f OD) 

MESSAG  (WHEN, 10,1  .  .'y.  3) 


ABUT 


ABUT') 


1  16 

lie 

115 


120 
130 
133 
150 


CALL  MFSSAG(TME0D,10,'AbUT','ABuT') 

CALL  EMOPL(G) 

CONTINUE 

CONTINUE 

CONTINUE 

IF (ILAST.NE  .1  )  60  TO  10 

CALL  OONEPL 

F0R>«AT(lX,IA,lX,F6.1,3(F6.2).6(F9,3).t(f7.3)) 

FOffWAT   (1X,I3»  4(F1P.1 ) ,7( F10.3)  ) 

FORf^AT  (1  X  ♦14F9.2) 

F0RI*AT  () 

END 
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s 

INTC 

D 

C 

•  0 

c 

*R 
«R 
*A 

I 
A 

E 
R 
R 
R 
A 
E 
R 
R 
T 
I 
I 
D 
N 
3-0 

ORI 

RT  ' 

ROOT 

R200 

I 

I 

I 

I 

I 

R 

VELO 

V 

T 
T 

T 
T 
R 
R 

X 
Y 
2 

X 

Y 
7 

V 
V 
V 
V 
V 
V 


UPROU 

ON  IN 

I^ENS 

OWMON 

ELT,N 

0»*WON 

0L4D6 

AN6TM 

NA  ,AN 

TCS  =  1 

2  =  0. 

L=0. 

0L1=0 

0L2-0 

OLA=0 

20EG  = 

LOEG  = 

OLIOG 

0L206 

1AN6  = 

6  =  0 

62=0 

ELTsD 

PRF  =  1 

TARGE 

GIN  A 

TRUE 
T  =  T 
TT  = 
F(IR, 
F(IR. 
F<IR. 
FdOS 
FCIOS 
200TT 
C  «  T 
ELOC  = 

TRAD 
RA01- 
RA02  = 

TRAT 
RAT1  = 
flAT2  = 
ANGT  = 
ANGTP4 

-R  = 
R  =  RT 
R=0, 
R=0. 

-T  s 
T  =  0. 
T  =  0. 
T  =  0. 

V-T 
XT  =  VE 
YT  =  0. 
2T  =  0. 
X  T  W  =  V 
YT^sV 
2T«  =  V 


TINE 
ITILI 
ION  P 

XT  ,Y 
PRF  ,T 
/INTC 

♦  T6S 
tXR  ,Y 
E»PAU 


INTCON 

2ES  TRAJECTORY  VARIABLES 

AUW(3) 

T,2T,V 

RAT2,I 

/     A2 

1  .T6S2 

ft.2R,V 

V 


XT,VYT,V2T,VEL0C,TlAN6,T2AN6,DELR,I6.TRATlt 
I,NPH»R90,1G2,ICTS.PI.IR,IH,I0S 

tELfROLl,ROL2,ROL^.A2DEG,FL0E6.ROLlD6  ,R0L206» 
,RT,R0OTT,R2d0TT,     DT,TRADlfTRAD?,RAN6T, 
XTM,VYT«,V2TM,WR20,AXT,AYT,A2T,AXT6,AYT6,A2T6, 


0. 
0. 
=  0. 
=  0. 
0. 


T  TRA 
T  INI 
RANG 
RUE  V 
TRUE 
EQ.1) 
EQ.2) 
EQ.3) 
•  EQ.1 
.EQ.2 
=  0. 
ARGET 
-ROOT 

-  =  T 
VELOC 
VELOC 

-  =  T 
T6S1* 
TGS2* 
RT 

=  RAN6 
RAOA 


JECTOR 

TIAL  T 

E  RELA 

ELOCIT 

ACCFLF 

RT  =  11 

RT  =  10 

RT  =  15 

)  ROOT 

)  ROOT 

VELOC 
T 

URN  RA 
•VELOC 
*VELOC 
URN  RA 
10.725 
10.725 

T-(RDO 
R  LOCA 


Y  DA 
AR6E 
TIVE 

Y  RE 
RATI 
275. 
00. 
0. 

T  =  -1 
T=-1 


TA 

T  POSITION 

TO  RADAR 
LATIVE  TO  RADAR 
ON  RELATIVE  TO  RADAR 


75. 

00. 


ITY  RELATIVE  TO  ORIGIN 

DIUS  OF  T6S-  6.  TURN 

/T6S1*10.7252 

/T6S2*10.7252 

TE    OF  TGS-  G.  TURN   (PAD/SEC) 

2/VELOC 

2/VELOC 

TT*OELT*FLOAT(NPRF)) 
TION  AT  T(C) 


TRUE  X,V,2  TARGET  POSITION 


=  TRUE  X,Y,2  VELOCITY 
LOC 


XT 
YT 
2T 


^3 


VR?D=R0OT7 

A-T  =  Tf<Uf 
AXT=0. 


X.YtZ  ACCELFRA7I0N 


=  TRUE  X,V,Z  ACCELERATION  (6'S) 


AyT=0. 

AZT=0. 
A-T& 

AXTG=0. 

A YT6=0. 

AZT6=0. 
ANA  =  VERTICAL  (AZ)  ACCELERATION  RELATIVE  TO  TARGET  ORIENTATION 
ANE  =  HORIZONTAL  (EL)  ACCELERATION  RELATIVE  TO  TARGET  ORIENTATION 

ANA^O, 

ANe=0. 
PAUV(I)  =  PITCH  AXIS  UNIT  VECTOR 
PAUVd)  =  X 
PAUV(?)  =  Z 
PAUV(3)  =  Y 

PAUV(1)=0. 

PAUV(2)=1 . 

PAUV(3)=0. 

END 
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FLLCS  ,S 


^''  V  R  ,  r< ,  P  V  k  P  ,  M  y/  R  T 


L  ,  U  ,  1    TO   L.. ,  -  A  ,  P  ) 
RP  r  ROLL  RATE  CF  AIRCRAFT  (DFG/StC) 
RPR  z    POLL  RATE  OF  AIRCRAFT  (PAD/SEC) 
ICTS  =  COORDINATED  TURN  SWITCH 

DTrDFLT 

RRrSE. 

PRRrPfi  *PI/ 180, 

CONDITIONAL 

( II .LT.NPH) 

INITIAL  TRAJECTORY 

ITCSrO 

XTIXT-»DELR 

•MEASURED'  AIRCRAFT  POLL,  PRECEEDIN6  MANUVFft 

ROL'*  =  ROL*»*  CkRR«DT  ) 
..FIN 
Tl  ANG.LT .R9Q ) 

SECOND   TRAJFCTORY 
R0L'*-P0L*4*  (  F'PR'i'DT  ) 

IF  (ROLM.GT.  ]  .55) 
.   R0LUrl.«^5 

...FIN 

IF  (  IG.GT.  i<^  ) 

I  C  T  S  —  "^ 

»ME*SUREr.'     flIPCRAFT     ROLL,     PPECFECING    MANUVER 

,        POL**  =  ROLH-lRRR*DT) 

..  .FIN 

IG-IG*! 

-     ANGLE     OF     TARGET 
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Hl^L  TKfiJtCTCr^Y 

t'LiSUfLr'     AlRCI^AFT     POLL 
F  (  h  0  L  «♦  .  G  T  .  C  .  C  ) 

R0Lt4  =  PCL4-(PRR*DT) 

..Fir, 

F  (ROL  ti.L  T  .  ..b  ) 

ROLUzr.o 

.  .  F  U.' 

T^Ar.&:TRATr«riLLT«FLCAT(NiFRF)«IG 

YT=YT-»r.'ELK*SlN(  T2ANG) 
ZTz:T*0EL(-*COS<T:ANe) 

VY  T=VFLOC*S  IN { T2A\6) 
V<.^  Trv£LOC*COS  (  T^A  NG^ 
FI\ 


.  .  F  T  N 
RETU^'N 
END 


CS    VERS10\    PJ2,6) 


qfi 


SUBROUTINE  FTkALI   < K  ,  X  .  m  ,  6  ,  Y  .  S  ,  Q  .  R  ,  P  .  I  N  ,  I  S  .  I  L  .  N , « 1  ,  l  ,  T 1  ,  7 2  , 
FTKALI S 


IT,T3,T)(  ,1ER) 


FUNCTION 
USA6F 

PARA«*ETERS 


IN 

IS 

IL 

N 

«1 

L 

T1 

T? 

IT 

T3 
TX 

lER 


ITERATE 
CALL  FT 

T1  ,12 
INPUT  S 

VECTO 

OF  TH 

CALCU 

X  AS 
INPUT/0 

X  IS 

OUTPU 

VECTO 
INPUT  m 

TRANS 
INPUT  M 
INPUT  0 
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